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ABSTRACT 


Aims. We study the evolution of stellar rotation and wind properties for low-mass main-sequence stars. Our aim is to 
use rotational evolution models to constrain the mass loss rates in stellar winds and to predict how their properties 
evolve with time on the main-sequence. 

Methods. We construct a rotational evolution model that is driven by observed rotational distributions of young stellar 
clusters. Fitting the free parameters in our model allows us to predict how wind mass loss rate depends on stellar mass, 
radius, and rotation. We couple the results to the wind model developed in Paper I of this series to predict how wind 
properties evolve on the main-sequence. 

Results. We estimate that wind mass loss rate scales with stellar parameters as M* oc R*fii ,33 MT 3 ' 36 . We estimate 
that at young ages, the solar wind likely had a mass loss rate that is an order of magnitude higher than that of the 
current solar wind. This leads to the wind having a higher density at younger ages; however, the magnitude of this 
change depends strongly on how we scale wind temperature. Due to the spread in rotation rates, young stars show a 
large range of wind properties at a given age. This spread in wind properties disappears as the stars age. 

Conclusions. There is a large uncertainty in our knowledge of the evolution of stellar winds on the main-sequence, 
due both to our lack of knowledge of stellar winds and the large spread in rotation rates at young ages. Given the 
sensitivity of planetary atmospheres to stellar wind and radiation conditions, these uncertainties can be significant for 
our understanding of the evolution of planetary environments. 


1. Introduction 

Low-mass main-sequence stars slowly lose mass and angu¬ 
lar momentum through magnetised stellar winds, though 
the properties of such winds are currently not well under¬ 
stood. These winds influence stellar rotation by removing 
angular momentum, causing their rotation t o slow down as 
they age ( Weber fe Davis|[l967 Kraft|1967 1 and leading to 
older stars having weaker magnetic fields and lower levels 


of high energy radiation (Giidel et al. 


2005 


Reiners 2012 Vidotto et al.|2014a l. The evolution of 


1997; Ribas et al. 


planetary atmospheres is highly sensitive to the surround¬ 
ing stellar environments. For example, high levels of X-ray 
and EUV radiation incident on a planetary atmosphere lead 


still remain (for su mm aries, see Krishnamurthi et al.||1997 
Herbst et al.| 2007 and Bouvier et al. ||2013|) . Rotation peri- 
ods are now known for tens of thousands of stars, with par¬ 
ticularly interesting studies being surveys of rotation peri¬ 
ods in stellar clusters of different ages. Stars start out their 
lives on the main-sequence with rotation rates that span 
two orders of magnitude, as can b e seen in young clusters 
such as NG C 2547 ©40 Myr ; Irwin et al. 20081 and Pleiades 
©120 Myr; Hartman et al. 2010| . This spread in rotation 
rates can be traced back to the early pre-main sequence and 
has the appearance that most stars lie on one of two district 


branches (Soderblom et al. 1993 Barnes 2003). As stars age, 
their rotation rates quickly converge, as can be seen in older 


Lammer et al. 

2003 

Tian et al. 

2008 


clusters such as M37 ©550 Myr; 

Hartman et al. 2009) and 

NGC 6811 ©1 Gyr;|Meibom et a 

. 2011). This convergence 

happens quicker for higher mass s 

;ars (Krishnamurthi et al. 


the same time, strong winds compress the magnetosphere 1997J) . At 100 Myr, most solar mass stars already lie on the 
and can cause non-thermal erosion of the atmosphere a sig¬ 
nificant amount of at mo spheric gas lies a b ove the magneto- 
spheric obstacle (Lichtenegger et al. , 20l0l |Kislyako va et al. 

|2014l. Therefore, the spin down of stars Is highly signifi- 
cant in the study of the evolution of planetary atmospheres 
and the develo pment of habitable planetary environments 
(Guedel et al. 2014). Stellar winds can therefore severely 
influence the atmospheres of planets, both directly and in¬ 
directly through their influence on the rotation rates of their 
host stars. 


slowly rotating track, and by 500 Myr, almost all solar mass 
stars have converged. On the other hand, for 0.5 M 0 stars, 
rotation rates are much more distributed at 100 Myr, and 
much less convergence has taken place by 500 Myr. Go¬ 
ing to masses below 0.35 M 0 , there is evidence that the 
distribution of rotation rates in fact broadens between the 
ZAMS and 10 Gyrs ( Irwin et al. [20 11 1, though this is diffi- 


Since the initial confirmation o f the s pin down of main- 
sequence stars as they age (Kraft 1967), a detailed obser¬ 
vational understanding of the rotational evolution of main- 
sequence stars has been achieved, though many questions 


cult to study due to the difficulty m measuring stellar ages 
for non-cluster stars. We do not consider such low stellar 
masses in this paper. After the convergence of stellar rota¬ 
tion rates there exists a one-to-one relation between stellar 


by © oc t 0 5 (Skumanich 

1972 

Barnes 2007 

Mamajck & 

Hillenbrand 2008). 
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A likely consequence of the spin down of stars and the 
decay of their magnetic dynamos is that their winds evolve 
in time. The link between magnetic activity and winds is 
poorly understood, so constraining the exact dependence of 
winds on stellar mass and rotation is difficult. Intuitively, 
we would expect that the solar wind mass loss rate was 
higher in the past as a result of stronger magnetic activ¬ 
ity. Observationally, there is some evidence for this (Wood 
et al.||2002 Wood et al.|2005 1, though the current observa- 
tional picture is unclear. There have been some theoretical 
studies of the past solar wind properties. |Belcher & Mac-| 
Gregor (1976) assumed that the processes responsible for 
the acceleration of the current solar wind, which they called 
simply ‘thermal’, have remained constant with age, and es¬ 
timated that the early solar wind could have had velocities 
of ~4000 km s -1 simply due to magneto-rotational acceler¬ 
ation. However, they found little influence on the mass loss 
rate. In reality, the ‘thermal’ acceleration is unlikely to be 
constant, but will depend on rotationally driven changes in 
the star’s magnetic field. More recent models for the evolu¬ 
tion of the solar wind have concentrated on the influences 
of these mech anisms flGriefimeier et al.||2004| Holzwarth & 
Jardine|2007 Sterenborg et al.|2011| Suzuki et al.|2013|). A 

’ els is thi 


general property of these models is that the mass flux in the 

solar wind was higher at previous times than it is currently, 
though it is unclear by how much this was the case. What 
is also unclear is how the wind speed has changed in time. 
Observations of the coronae of young solar analogues show 
that coronal temperatures decay with magnetic activity as 


stars spin down (Giidel et al. 1997 Telleschi et al. 


2005), 


which suggests that the solar wind temperatures, and there- 
fore speeds, were higher at young ages. However, the link 
between coronal and wind properties is currently unclear, 
and so we should be cautious when making such conclu¬ 
sions (see Section 4.1 of Paper I of this series for a detailed 
discussion on this issue). 

Another open question about stellar winds is how they 
depend on stellar mass. Naively, we would expect that lower 
mass stars have lower mass loss rates due simply to the 
smaller surface areas. However, the situation is unlikely to 
be so simple; for example, at a given rotation rate, low- 
mass stars have higher X-ray surface flu xes, f’x■ th an high- 
mass stars in the unsaturated regime. Vidotto et al. (2011) 
used 3D magnetohydrodynamic (MHD) simulations of the 
wind from the ~0.3 M 0 rapid-rotator V374 Peg and argued 
that the star is li kely to have a more powerful wind than 
the current Sun. Cohen et al. (2014) modelled the wind 
of the ^0.3 M 0 moderate-rotator EV Lac and derived a 
mass loss rate per unit surface area that is an order of 
magnitude above the current solar wind. A significant open 
question when it comes to the mass dependence of stellar 
winds is the possible saturation of mass loss rates at fast 
rotation. In X-rays, the saturation threshold for low-mass 
stars is lower than for high-mass stars, meaning that the 
most rapidly rotating low-mass stars are not able to become 
as active as the most rapidly rotating high-mass stars. If 
something similar is happening for winds, then the mass 
loss rate could have significantly different mass dependences 
in the saturated and unsaturated regimes. 

In this paper, we study the rotational evolution of low- 
mass main-sequence stars and the influence that this has 
on the properties of stellar winds. We concentrate on rota¬ 
tional evolution for two reasons: firstly, by fitting our ro¬ 
tational evolution model to the observational constraints, 


we attempt to derive a scaling law for wind mass loss rates 
as a function of stellar parameters, and secondly, we use 
rotational evolution to predict the time evolution of wind 
properties on the main-sequence. This study is made pos¬ 
sible by the recent derivation of the dependence of wind 


torque on mass loss rate (Matt et al. 20121, and a recent 


tion in young stellar clusters (e.g. Irwin et al. 

2007 

Hart- 

man et al. 2009| |Irwin et al. 2009 

|Meibom et al. 

2009; 

Hartman et al. 

2010). In Section |2 

we develop our rota- 


porated into the model as a free parameter. In Section [3j 
we collect observed rotation periods for over 2000 stars in 
seven young clusters to derive observational constraints on 
the rotational evolution of stars on the main-sequence. In 
Section [4j we fit the free parameters in our rotational evo¬ 
lution model to the observational constraints, allowing us 
to derive wind mass loss rates as a function of stellar mass, 
radius, and rotation rate. In Section [5] we summarise the 
results of our rotational evolution model. 

The results of our rotational evolution models are also 
that we can predict both the rotational evolution of individ¬ 
ual stars and the distributions of rotation rates at the ages 
and stellar masses that we consider. In the second part of 
this paper, we couple the results of the rotational evolution 
model to our wind model derived in Paper I. In Section [6j 
we show how stellar wind properties, such as mass fluxes 
and wind speeds, evolve on the main-sequence. Finally, in 
Section [TJ we summarise our results and discuss their im¬ 
plications for the evolution of planetary habitability. 


2. Rotational Evolution: the Physical Model 

In order to study the evolution of stellar winds from stars 
of different masses and ZAMS rotation rates, we first need 
a good understanding of the rotational evolution of main- 
sequence stars. For this purpose, we develop a rotational 
evolution model for stars of masses between 0.4 M 0 and 
1.1 M 0 and between ages of 100 Myr and 5 Gyrs. In order 
to constrain the free parameters in our rotational evolu¬ 
tion model, we collect measured rotation periods for several 
young clusters from the literature. 


As summarised by Krishnamurthi et al. (1997), the basic 
ingredients of any rotational evolution model are the initial 
rotation rate, the internal structure of the star, and the rate 
at which angular momentum is removed from the star by 
the wind. Since we do not consider rotational evolution at 
ages younger than 100 Myr, we constrain our initial condi¬ 
tions using observations of rotation in young clusters with 
ages close to 100 Myr. To determine the internal structure 
of the star, and most importantly, the moments of inertia as 
a function of stellar mass and age, /*(M*, t ), we run stellar 
evolution models for each stellar mass of interest. For the 
stellar evolution calculations, we use the ID hydrodynamic 
stella r evolut ion c ode described in|Heg er et al. (20001, Brott 
et al. (2011aI, and Brott et al. (2011b). For the simulations, 
we assume solar metallicity and neglect the effects of stel¬ 
lar rotation and mass loss (which would anyway only have 
negligible effects on the models). In addition to the stellar 
moments of inertia, we use these simulations to derive stel¬ 
lar radii and bolometric luminosities as a function of stellar 
mass and age, which we use as important input in the rest 
of the calculations in this paper. 
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An important ingredient in several rotational evolution 
models is core-envelope decoupling. Core-envelope decou¬ 
pling happens when the timescale for angular momentum 
transport within the star is not negligible compared to the 
timescale over which the moment of inertia of the star 
changes and the timescale over which angular momentum 
is removed from the stellar surface by the wind. The result 
would be that stars arrive at the ZAMS with cores that ro¬ 
tate more rapidly than the surfaces. Several previous stud¬ 
ies have included a core-envelope coupling timescale as a 
free parameter in their models and found that short cou¬ 
pling timescales of a few 10 Myr are required to reproduce 


the observational constrai nts (e.g. MacGregor & Brenner 
1991 Spada et al. 2011 Gallet & Bouvier 2013[), which 


leads to moderate differential rotation between the core and 
the envelope. In order to properly take into account core¬ 
envelope decoupling, we would ideally have to model rota¬ 
tional evolution before the ZAMS in order to determine the 
initial core rotation rates for ZAMS stars. Given the short 
coupling timescales found in these previous studies and the 
fact that we are able to fit our models to the observational 
constraints without assuming core-envelope decoupling, for 


where the wind torque has a mass loss rate dependence), 
and the internal structure of the star. In order to explain 
the fast spin down of young rapidly rotating stars, it is 
necessary that their winds carry awaj_ angular momentum 
faster than the current solar wind (iBelcher Sz MacGregor 
1976), which can easily be explained by higher magnetic 
fiem strengths and mass loss rates. In general, Skumanich 
style spin down requires that the wind torque depend on 

However, it is difficult to explain the slow spin down of 
the most rapidly rotating stars with such a strong depen¬ 
dence. Instead, a much weaker dependence between wind 
torq ue and rotation above a certa i n threshold is required 
(e.g. MacGregor fe Bren ner 1991). Barnes & Sofia (19961 
and Bouvier et al. (|1997|) showed that this threshold must 
be lower for "low-mass stars in order to explain the slower 
spin down of rapidly rotating low-mass stars. The reason for 
this threshold is likely the saturations of the wind mass loss 
rates and magnetic field strengths at high rotation rates. 

In the absence of core-envelope decoupling, the rate at 
which a star spins down is given by 


lope 


simplicity, we assume solid-body rotatio: 

The fact that a rotating magnetised star that loses mass 
through an ionised winds will spin down was recognised 
immediately af ter the first predictions of the existence of 
(Parker||1958 |Schatzman||1962 ). 
lown is that the wind has a lai 


dt 


1 ( d, ‘n 


(1) 


the solar wind 


The reason 
larger specific 


for this spin down 
angular momentum than the material in the star, mostly 
due to the angular momentum contained i n the stresses 
of the magnetic field (Weber & Davis 1967|). As the wind 


propagates outwards, the angular momentum held in the 


where /* is the star’s moment of inertia and r w = dJ/dt is 
the torque on the star by the wind, where J is the star’s 
angular momentum. On the pre-main-sequence, the term 
involving dl+/dt is very important; on the main-sequence, 
this term is negligible. Therefore, in order to predict how 
a star’s rotation evolves in time, we require a model for 
calculating the wind torque. 


magnetic field is then transferred to the gas (see Fig. 3 of 

Matt et al. 

( 

Vidotto et al.|2014b 

). Another way to think about it is that 

lar mass, M 



where the torque is related to the 


per unit mass of the wind in the equatorial plane at the 
surface of the star is n±R+. However, in the presence of 
a magnetic field, the wind torque is equivalent to what it 
would be if the material was held in strict corotation with 
the star out to the Alfven radius, Ra, and then released, 
and therefore the angular momentum per unit mass lost in 
the wind in the equatorial plane is tt+R^. 

The rate at which stars lose angular momentum depends 
on the stellar magnetic field, the wind mass loss rate, the 
stellar mass and radius, and the angular velocity of the star. 
Relating the wind torque to these parameters is difficult. 
Most rotational evolution models have used a formula for 


fl*, magnetic field strength, B^ip, and mass loss rate, M „, 


as 


t' = KfB^Ml- 2 m Rf n+2 




(K%vi c + niRiy 


( 2 ) 


calculating wind torques derived by Kawaler (1988) based 
on simple ID approximations for the wind an a the magnetic 
field. In most formulations of the Kawaler (1988) formula 
that have been used, the wind torque does not depend on 
the mass loss rate, which is not reasonable. A mo re realistic 
formula for wind torques was derived by Matt et al. (2012) 
using a grid of 2D MHD wind models. The main advantage 
of this approach is that simplifying assumptions about the 
Alfven surface do not need to be made since this is solved 
in the MHD simulations self-consistently. 

Formulating rotational evolution models using these 
torque formulae is a tricky business that requires detailed 
knowledge of the field strengths, mass loss rates (in models 


1 Recently, Gallet & Bouvier (20151 argued that low mass stars 
(«0.5 M 0 ) have much longer core-envelope coupling timescales 
of >100 Myr. 


where I\\ = 1.3, K 2 = 0.0506, m = 0.2177, v esc is the 
surface escape velocity, and all quantities should be in cgs 
units. This torque formula is calculated from the results of 
50 MHD wind simulations and implies that when all other 
parameters are held constant, the wind torque depends on 
M* " 56 . The fact that the M* dependence should be weaker 
than linear can be seen from simpler considerations. The 
torque associated with the wind in the equatorial plane is 
given by t w = Af*fThe mass flux in a wind can be in¬ 
creased either by increasing the wind speed or by increasing 
the wind density. Increasing the mass loss rate while keep¬ 
ing the magnetic field strength constant leads to a lower 
value of the Alfven radius, either because the wind acceler¬ 
ates faster or because the higher wind density leads to lower 
Alfven speeds. Therefore, increasing the mass loss rate leads 
to a slower than linear increase in the wind torque. 

This formulation represents a significant advance in our 
ability to predict stellar wind torques and is li kely more re¬ 
alistic than the torque formula from Kawaler (1988). How¬ 
ever, as we discuss in more de tail in Section |4.3[ the sim¬ 
ulations of Matt et al. (2012) did not take into account 
all of the detailed physical mechanisms operating in stellar 
winds. Furthermore, our constraints on the values of Hdip 
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and M* are likely to contain uncertainties. To take into ac¬ 
count these uncertainties, we add an extra variable into the 
torque formula, such that the actual torque that we use in 
our model, r w , is related to the torques calculated using 
Eqn. [2] by 


T'w bh T . 


(3) 


where K r is a free parameter in our model. A similar fac- 
tor was introduced into Eqn. 0 of |Gallet fe Bouvier ( )2013| ) 
(though they did so by changing the value of K\), who 
found I\ T « 2. 

In order to calculate the torque on a star from the wind, 
we therefore need to be able to calculate M* and I?dip- We 
derive the mass loss rates by fitting our rotational evolu¬ 
tion model to the observational constraints. We make the 
assumption that in the unsaturated regime, the mass loss 
rate per unit surface area of a star has power-law depen¬ 
dences on its mass and rotation rate, such that 


—jimnr 


(4) 


201 4). A measure of magnetic field strength that is likely 
to be more reliable is X-ray emission. Combining the rela- 


of Lx oc $k 15 from 

Pevtsov et al. 

(2003 

with the relation 

between Kx (where Kx = Lx/^bo 

) and Rossby number of 

Rx oc P ro 218 from 

Wright et al. 

(2011 

) would imply for 

a given stellar mass that fB oc J 

p—i.y 
rot, ? 

though using the 


where a and b are free parameters. This assumption is likely 
a simplification; in real stellar winds it might not be possible 
to genuinely separate the stellar mass and rotation depen¬ 
dences in this way. We discuss our fits to these parameters 
in detail in Section [4j We choose to include a surface area 
dependence in the above assumption for mass loss rate since 
this is likely to be more realistic. However, this is unlikely to 
influence our results significantly given that we determine 
the strength of the mass dependence of M* by fitting the 
rotational evolution models. Removing the surface area de¬ 
pendence would simply lead to us deriving a different value 
of b. As discussed below, we assume that the dependence of 
M* on 12* saturates at fast rotation. 

How the magnetic field strength should be calculat ed is 
a difficult issue. It has been known since Kraft (1967) and 


steeper relation of Rx oc P rot from Wright et al. (20111 
would imply that fB oc P“ 2 ' 3 . 

We stress, however, that the magnetic field term in the 
torque formula (Eqn. pi is the equatorial strength of the 
dipole component of the field and does not correspond to 
the surface average field strength. Naively, we would expect 
that the strengths of the different components of the field 
scale with rotation in the same way. However, using the re¬ 
sults of recen t Zeem a n-Dop pler Imaging studies of over 70 
stars, Vidotto et al. (2014a) showed that the large scale field 
strength averaged over the stellar surfac^J By, scales with 
rotation as By oc P“ t l:32 , which is weaker than the relations 
found for fB. Although this is not the dipole component of 
the field, there is an approximately linear scaling between 
By and Pdip in the sample of stars (Aline Vidotto, private 
communication). One thing that is unclear is whether the 
magnetic field strength scales with rotation rate directly or 
with some other parameter, such as Rossby number. We 
consider the correlation betwee n By and Rossby number 
shown in Vidotto et al. (2014a I to be better than the cor¬ 
relation between By and rotation period (compare their 
Fig. 3 and Fig. 4), and so we assume that Pdip scales with 
Rossby number. Therefore, when 12* < 12 sat , we get 


Pdip = B , 


dip,© 


/ 12* r* \ 

K^QTq) 


1.32 


(5) 


Skumanich (11972) that rapidly rotating stars~have higher 
levels of magnetic activity than slowly rotating stars. At 
a given stellar mass, the magnetic field strength has ap¬ 
proximately a power law dependence on rotation rate, such 
that P* oc 12* fo r slow rotators, and saturates at fast ro¬ 
tation (e.g. Reiners|2012 1. However, exact quantification of 
this relation is difficult and requires knowledge of the index 
in the power law, the constant of proportionality, and the 
(mass-dependent) location of the saturation threshold. 

Most studies of magnetic field strengths in the past have 
focused on the surface average field strength, often denoted 
fB, where / is the surface magnetic fil ling factor and B 
is the magnetic field strength. Saar (1996) and Saar (2001) 
found that fB is proportional to P“("' and Po~ 1,2 respec¬ 
tively. The Rossby number, Ro, is given by Ro = Prot/t*, 
where r* is the convective turnover time. In the latter case 
especially, this relation is likely to be made artificially shal¬ 
lower due to the inclusion of stars in the saturated regime 
and the lack of real measurements for slow rotators (the 
measurements for two of the three stars slo wer than 12 days 
were upper limits in the sample of Saar| 1996). At the same 
time, measurements of fB based on Zeeman broadening 
are very difficult for weak magnetic fields and can be dif- 
ficult to distinguish from rotational broadening (Reiners 


and when f2* > f2 sat , the dipole field strength remains at 
the saturation value. The convective turnover times, r* and 
70, should not be confused with the wind torque, r w . 

A very important part of rotational evolution models is 
the saturation of the magnetic field strength and the mass 
loss rate. In the saturated regime, the dependence of the 
wind torque on 12* must be much weaker than in the unsatu¬ 
rated regime. Without saturation, the most rapidly rotating 
stars would spin down much faster than is observed. Based 
on the slow rotational evolution of young rapidly rotating 
stars, there can be no serious doubt that both mass loss rate 
and the dipole component of the field saturate at fast rota¬ 
tion. However, it is not clear that saturation takes place at 
the exact same rotation rates for the mass loss rate and the 
magnetic field. Similarly, it is not clear that saturation takes 
place at the same rotation rates for different components of 
the magnetic field (e.g. Pdipi fB, etc.), and for indicators 
of magnetic activity such as X-ray emission. For simplicity, 
we assume that both M* and P^ip saturate at the same ro¬ 
tation rate for a given stellar mass. The exact rotation rate 
at which saturation takes place is not well constrained, al¬ 
though good hints on the saturati on level can be derived 
from X-ray emission. For example, Pizzolato et al.| (20031 
estimated that the rotation period for saturation depends 
on stellar luminosity as P sat ~ 1.2 (Lboi/P©) -1 ^ 2 ) where 


2 The subscript V here represents magnetic field strengths that 
have been derived from Stokes V measurements, as opposed to 
Stokes I measurements which are used to derive fB. 
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P sat is in days. A similar result of P sat « 1.6 (Pboi/ T 0 )~ 1//2 
was estimate d using a l arger sample of stars by |Reiners| 
et al. (20141. Similarly, Wright et al. ( 2011| ) determined 


that saturation in X-ray emission takes place for all stars 
at approximately a Rossby number of 0.13. How the satu¬ 
ration rotation rate is then determined from the saturation 
Rossby number depends on the exact form of the convec¬ 
tive turn over ti m e. The convective turnover times derived 
by Wright et al. (20111 imply that for 1 M 0 and 0.5 M 0 
stars, saturation occurs at rotation periods of 1.9 days and 
4.7 days respectively. In our model, we make the assump¬ 
tion that the saturation rotation rate for the magnetic held 
and the mass loss rate has a power law dependence on stel¬ 
lar mass, 


flsat (HP) — H sa t 50 



( 6 ) 


where f l 


sat,© 


is the saturation rotation rate for solar mass 


stars . Fro m the r esults of Pizzolato et al. (2003) and Rein- 
(20141, we would expect that sT 


ers et al. 
and 17 SR 


^sat,© 


21 Hr. 


respectively. On the other hand, the results of 


Wright et al. (2011) suggest that H sat , 0 is 14.5 f2 0 . In our 
models, we make the assumption that fl sat 0 is 15 O 0 and 


treat c as a free parameter. Th e results of|Pizzolato et ah] 


(2003) and 

Reiners et al. 

)2014 

i suggest that c — 2. On the 

other hand, the results of 

Wrig 

ht et al. 

(2011) suggest that 


parameters in our model by fitting them to observational 
constraints discussed in Section [Sj We find that in order to 
properly reproduce the rotational evolution of lower mass 
stars, values of c larger than 2 are required, and we there¬ 
fore set c = 2.3 in our rotational evolution model. We use 
this value since it is the closest value to the above X-ray 
constraints that allows us to get visually acceptable fits to 
the observational rotational evolution. 


Clearly, to use Eqn. [5] we need to know the equatorial 
strength of the dipole component of the Sun’s magnetic 
field. This is difficult because the dipole component on the 
Sun varies over the solar cycle, with a maximum strength of 
a few G around cycle minimum. Therefore, given that the 
mass loss rate of the solar wind is approximately constant 
over the solar cycle, we expect that the wind torque on the 
Sun has a time dependence due to changes in the dipole 
field strength. To see how the strength of the dipole com¬ 
ponent of the field varies over the solar cycle, we use dipole 
field strengths derived from magnetic maps of the solar sur¬ 
face by the Wilcox Solar Observatory (WSO)[^] WSO has 
been measuring the solar magnetic field daily for almost 
forty years, and has derived synoptic magnetic maps for 
Cycle 21, Cycle 22, and Cycle 23. In Fig. [T] we show their 
measurements of the dipole held strength over these three 
cycles. Since we are interested in the evolution of stellar ro¬ 
tation rates over hundreds of millions of years, we are only 
interested in the average torque exerted on the star, and not 
variations over short time scales. Since the torque is pro¬ 
portional to we calculate the mean value of B°jp 7 Q 

over these cycles, which corresponds to i?dip,o ~ 1-35 G, 
which we then take as i?dip,© in Eqn. pi Since the R°jp 7 0 
term acts simply as a multiplicative constant in Eqn. [2] the 
exact value of i?dip,© does not matter and changing it would 
just lead to a different value of the free parameter K T . 



Fig. 1. Plot showing the time variations in the equatorial 
strength of the dipole component of the solar magnetic field 
as measured by the Wilcox Solar Observatory. The dashed hor¬ 
izontal line shows Bdip = 1-35 G. The dotted vertical lines show 
the beginnings of the previous three activity cycles. 


We now have a rotational evolution model with four free 
parameters. These are K T in Eqn. [3j a and b in Eqn. [4] and 
c in Eqn. [6] We constrain these parameters in Section [4] 


3. Rotational Evolution: Observational Constraints 


To constrain the free parameters in the rotational evolu¬ 
tion model, and to study the distribution of stellar ro¬ 
tation rates (and therefore wind properties) at any given 
age, we collect measured rotation periods from the litera¬ 
ture. The data that we collect are for seven young clusters: 
these are the Pleiades (—125 Myr), M50 (—130 Myr), M35 
(—150 Myr), NGC 2516 (-150 Myr), M37 (-550 Myr), 
Praesepe (—580 Myr), and NGC 6811 (—1000 Myr). Al¬ 
though rotation data is available for other young clusters, 
we choose to use only these clusters due to the large num¬ 
ber of measured rotation periods available for these clusters. 
The only exception is NGC 6811, for which many fewer ro¬ 
tation periods are known, and all of which are for stars with 
masses near 1 M 0 . However, given the advanced age of this 
cluster, including it gives us essential information about the 
rotational distr ibution at 1 Gyr. As can clearly be seen in 
Fig. 1 of Bouvier et al. (2013), the rotational distributions 
in the other young clusters that have been studied are con¬ 
sistent with the rotational distributions that we consider. 
A significant uncertainty in our analysis comes from uncer¬ 
tainties in the ages of the clusters that we use. Determining 


the a 

ges of 

et al. 

2013) 


the given ages as correct. 

We follow the method of several previous stud ies (|Irwin 
et al.||2007l |Bouvier||2008l |Hartman et al.|[2009[ |Ga 


let & 


Bouvier||2013 1 by using percentiles of the rotational distri¬ 
butions to fit our rotational evolution models. As described 
below, we combine the measured rotation periods from stars 
in several clusters to get rotation period distributions at 


3 The data from the Wilcox Solar Observatory can be obtained 
from wso.stanford.edu 
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Fig. 2. Plots showing the distributions of stellar rotation rates at ~150 Myr (left panel) and at ~550 Myr (right panel). The 
distribution at ~150 Myr is a combination of measurements in the Pleiades, M50, M35, and NGC 2516 and the distribution at 
~550 Myr is a combination of measurements from M37 and Praesepe, as discussed in the text. Histograms showing the distributions 
of rotation periods in different mass bins for these two sets of stars are shown in Fig. [3] 




log(P rot ) (days) log(PJ (days) 

Fig. 3. Histograms showing distributions of stellar rotation rates at ~150 Myr (left panel) and at ~550 Myr (right panel) at 
three different mass bins. The stellar samples from which these distributions are derived are shown in Fig. [2] 


~150 Myr and ~550 Myr, and we use rotation data from a 
cluster at 1000 Myr to constrain the rotation of solar mass 
stars at this age. We bin the stars in these distributions by 
mass and then calculate the 10th, 50th, and 90th percentiles 
of the distribution of angular velocities, fl*, for each mass 
at each age. We produce bins for stellar masses of 0.50, 
0.75 and 1.00 M 0 , assuming a bin width of 0.2 M 0 (such 
that the 0.5 M 0 bin contains all stars with masses between 
0.40 M 0 and 0.60 M 0 ). The result is that for each stellar 
mass, we have time sequences of the rotation rates for stars 
at the 10th, 50th, and 90th percentiles of the rotational dis¬ 
tribution. We make the assumption that a star starting at 
a certain position in the rotational distribution remains at 
this position as the distribution evolves, which means the 
time series of a given percentile can be used as observa¬ 
tional constraints on the real spin down tracks of stars at 
that percentile. This is likely to be reasonable given that 
a star’s magnetic activity, and therefore the rate at which 
it spins down, is primarily determined by its rotation rate, 
which means that a star at a given rotation rate is unlikely 
to overtake a slower rotating star as it spins down. 
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The percentiles for each mass bin and each cluster are 
given in Table. [T] The uncertainties in each value are calcu¬ 
lated using the bootstrap method and are mostly very small 
given the large number of stars in each mass bin. However, 
these uncertainties should be taken with caution; they rep¬ 
resent the uncertainties in the percentiles due to the fact 
that they are calculated from a finite number of stars, and 
do not take into account errors from other sources such as 
uncertainties in the rotation period measurements, biases 
in the distributions of period measurements, uncertainties 
in the stellar mass determinations, and uncertainties in to 
what extent the time sequences of percentiles can truly be 
used as a proxy for the rotation histories of stars. We do 
not consider the uncertainties in the percentiles when cal¬ 
culating the best fit values for the free parameters in the 
rotational evolution model. 


The two youngest clusters in our sample are the Pleiades 
and M50 wit h approximate ages of 125 Myr and 1 3 0 Myr 
respectively ( Stauffer et al.|1998 Kalirai et al.|2003 ). Hart 


man et al. (2010|) derived rotation periods for 383 Pleiades 


stars, of which 332 are in the mass range of interest to us. 
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Mass Bin (Mq) 

n* 

D*(Dq) 

om (fi 0 ) 

^10 (fi 0 ) 

^50 (^ 0 ) 

n 90 (d 0 ) 

Pleiades + M50 + M35 + NGC 2516 







(~150 Myr) 







All 

1556 

26.00 ± 1.08 

42.73 ± 2.79 

3.13 ±0.06 

7.66 ±0.38 

74.01 ± 2.54 

0.50 

463 

30.62 ± 1.54 

33.14 ± 1.82 

3.61 ±0.27 

17.35 ±1.42 

80.58 ±3.37 

0.75 

573 

26.75 ±2.05 

48.92 ± 5.09 

2.97 ±0.07 

5.30 ±0.33 

75.38 ± 5.32 

1.00 

230 

20.27 ±2.87 

43.65 ± 6.97 

3.72 ±0.36 

6.66 ±0.26 

48.64 ±6.83 

M37 + Praesepe 
(~550 Myr) 

All 

639 

7.17 ±0.39 

10.02 ±0.67 

1.97 ±0.07 

3.49 ±0.07 

19.02 ± 1.86 

0.50 

119 

15.67 ± 1.52 

16.87 ± 1.17 

1.57 ±0.04 

8.41 ±3.00 

43.88 ± 5.05 

0.75 

249 

4.59 ±0.33 

5.17 ±0.75 

2.22 ±0.03 

3.00 ±0.07 

8.33 ± 1.41 

1.00 

134 

4.30 ±0.15 

1.74 ±0.47 

3.14 ±0.12 

3.96 ±0.09 

5.42 ± 0.69 

NGC 6811 
(~1000 Myr) 

All 

51 

2.69 ±0.12 

1.17 ±0.56 

2.41 ± 0.06 

2.53 ±0.02 

2.75 ±0.10 

1.00 

36 

2.75 ±0.16 

1.39 ±0.68 

2.31 ±0.08 

2.53 ±0.03 

2.83 ±0.96 


Table 1. Table showing statistical properties of the distributions of rotation rates that we use to constrain our rotational evolution 
model. All of the angular velocities are given in units of the solar angular velocity, fi©, which we define as the Carrington rotation 
rate of 2.67 x 10~ b rad s -1 . As described in the text, the first two distributions are constructed by combining measured rotation 
periods from different clusters with approximately identical ages. Both of these aggregate clusters are shown in Fig. [2] For each 
cluster, we bin the stars by mass assuming a bin width of 0.2 M© and calculate for each mass bin, the number of stars in each bin, 
n*, the mean angular velocity, H*, the standard deviation in the angular velocity, an, and the 10th, 50th, and 90th percentiles of 
the distributions, given by fiio, SJso, and flgo respectively. The values marked as ‘All’ are calculated from all stars with masses 
between 0.4 M@ and 1.1 M©. The errors on the percentiles are calculated using the bootstrap method. We do not report statistics 
for NGC 6811 below the highest mass bin because very few stars lie below 0.9 Mg. 


Rotation in M50 was studied by Irwin et al. (2009) who 
measured rotation periods for 812 stars, of which 601 stars 
are in the mass range of interest to us. At approximately 


150 Myr (von Hippel et al. 2002 

Jeffries et al. 

2001 

). Mei- 

bom et al. (2009| determined rotation periods 

or 441 stars 


to us and have reported (B — V)q values that we use to 
calculate the stell ar mass. The stars in the sample of |Mei-| 


bom et al. (2009) all have masses above 0.6 M 0 . Luckily, 
the measurements of M35 are nicely complemented by the 
r otati on per iods derived for stars in NGC 2516 by |Irwin| 


et al. (2010) provided a comparison of these four clusters 
and snowed that they have almost identical rotation period 
distributions (see Fig. 14 and Fig. 15 of their paper). This 
provides good support for the assumption that a time series 
of percentiles in the rotational distributions can be used as 
a proxy for the evolution of individual stars. By combining 
these clusters, we have a distribution of ~1500 stars which 
we assume represents rotation at 150 Myr. 

In Fig. [2} we show the distribution of rotation rates 
against stellar mass for stars at 150 Myr based on combining 
these four clusters. The structure of this distribution is well 


known in t he literature (e.g., Herbst et al. 2007 Bouvier 
et al. 20131. The most prominent feature is a track of slowly 


rotating stars, with higher mass stars on the track having 
faster rotation rates than lower mass stars. At 150 Myr, the 
track dominates for stars with masses above 0.8 M 0 and 
clearly extends down to stellar masses of 0.6 Mq. At later 
ages, this track dominates even more as more quickly ro¬ 


tating stars spin down (Bar nes| ;2003). Eventually the track 
contains almost all stars (excluding tight binaries that are 
tidally locked and possibly stars below 0.35 Mq). Another 
feature that is clearly evident in Fig. [2] is the track of rapid 
rotators, with a large distinct gap between the two tracks. 
The existence of the track of rapid rotators and the gap 
was seen in the Pleiades by van Leeuwen et al. (19871 and 
has been extensively discussed in the literature. The exis¬ 
tence of these two tracks can be traced back to the pre-main 
sequence. For example, two distinct tracks are clearly visi¬ 
ble at all masses in the distribution of rotation rates in the 


et al. (2007). They derived rotation periods for 362 stars, 
of which 140 are in the mass range of interest to us, all of 
which are below 0.7 Mq. It is quite remarkable that the 
distributions of rotation rates at every stellar mass derived 
for the se fo ur clusters are completely consistent. [Hartman 


~15 Myr old cluster h Per (Moraux et al. 2013). Histograms 
showing the distributions oFrotation periods in three dif¬ 
ferent mass bins are shown in Fig. [3] 

The next cluster in the age sequence that we u se is M37, 
at an age of ~550 Myr (Hartman et al. 20081. Hartman 


et al. (20091 reported rotation periods for 575 stars in M37, 
of which 495 stars are in the mass range of interest to us. 
Given its age and the large number of measured rotation 
periods, this cluster is probably the most important one 
for our study since it gives us reliable constraints on how 
quickly the rapidly rotating stars of different masses spin 
down. At approxim ately the s am e age is Praesepe, with an 
age of ^570 Myr (Delorme et al. 2011|). Rotation periods for 
stars in Praesepe w ere reported by~ De lorme et al 7| (2011), 
Scholz et al. (20111, and Agiieros et aT (|201ip, of which 108 
are in the mass range of interest to us. For our distribution 
of rotation rates at ~550 Myr, we combine the rotation 
periods collected for M37 and Praesepe, and assume an age 
of 550 Myr for the resulting set of 603 stars. The rotational 
distribution for stars in this aggregate cluster is shown in 
Fig- [2} 

For solar mass stars, the slowly rotating track seen at 
150 Myr has clearly spun down significantly and most of 
the rapid rotators have already converged onto the slowly 
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rotating track. The 10th percentile for the 1.0 M 0 mass 
bin decreases from 4 f2 0 to 3 fi Q between these two ages 
and the 90th percentile decreases from 50 fl 0 to 5 f2 0 . 
At lower masses, although the rotational distribution has 
clearly evolved towards slower rotation, the rapidly rotat¬ 
ing track seen at 150 Myr is still visible, and much less 
convergence of the two tracks has taken place. The 10th 
percentile for the 0.5 M 0 mass bin decreases from 4 fl 0 
to 2 f2 0 between these two ages and the 90th percentile 
decreases from 86 f2 0 to 44 fi 0 . Interestingly, the hole at 
intermediate masses and rotation rates seen at 150 Myr 
contains a large number of stars by 550 Myr. 

The final cluster that we consider is NGC 6811, at an age 
of ^1000 Myr. Rotation periods for 71 stars in NGC 6811 
were reported by Meibom et al. (2011), of which 51 stars 
are in the mass range of interest to us. Unfortunately the 
sample of stars does not contain any stars below 0.8 M 0 
and so cannot help us constrain the spin down timescales 
for the rapidly rotating low-mass stars. All of the stars have 
converged onto the slowly rotating track, and no rapid ro¬ 
tators exist. There is almost no difference between the 10th 
and 90th percentiles. 

Unfortunately, due to the difficulty of measuring the 
ages of individual stars, much less is known observation- 
ally about the evolution of rotation beyond 1 Gyr. Most of 
what is known is based on the current rotation rate and 
age of the Sun. It is often assumed in the literature that 
once all stars evolve onto the slowly rotating track seen in 
young clusters, this track spins down in a way that pre¬ 
serves its shape on the rotation-mass diagram, and there¬ 
fore, the rotation period as a function of stellar mass and 
age can be represented as the product of two functions, 


i.e. fiot( AAU) = f(M+)g(t) (Barnes 2003 Barnes 2007). 


Recently, Mamajek & Hillenbrand ( 2008| ) d erived an ex¬ 
pression for the shape of this track as a function of stellar 
(B — V) 0 and age, such that 


P rot = 0.407 [(B - V) 0 - 0.495] 0 ’ 325 f°' 566 , (7) 

where t is in Myr, P rot is in days, and (P — U) 0 represents 
the mass dependence. The shape of the mass dependence 
is based on the shapes of the slowly rotating tracks in the 
Pleiades (~125 Myr) and Hyades (~625 Myr), and the de¬ 
pendence on time is influenced strongly by the current ro¬ 
tation rate of the Sun. In order to constrain our rotational 
evolution models at older ages, we assume the above equa¬ 
tion is approximately true for stars beyond 1 Gyr. Using the 
above formula, we calculate the rotation period for stars in 
each of the six mass bins that we consider in our rotational 
evolution model at 5 Gyr. However, we warn that the shape 
of the mass dependence is not well constrained observation- 
ally beyond the age of Hyades. This represents a significant 
uncertainty in the analysis of the next section. 





4. Rotational Evolution: Determination of the Free 
Parameters 


In the previous sections, we construct a rotational evolu¬ 
tion model and collect rotation periods for more than 2000 
stars in several young clusters. In this section, we use the 
observational results to constrain the free parameters in the 
rotational evolution model. We evolve a star’s initial rota¬ 
tion rate forward in time by integrating Eqn. |T] using the 
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Fig. 4. Plots showing how our goodness-of-fit parameter, A', 
depends on the free parameters a and b. The upper panel shows 
X against a considering only the 1 M 0 mass bin. The vertical 
dashed line shows our best fit value of a ss 1.33. The middle 
panel shows X against b assuming the above best fit value for a 
and taking into account all three mass bins. The vertical dashed 
line shows our best fit value of 6 « —3.36. The lower panel shows 
a contour plot of X against both a and b, taking into account 
all mass bins. The red cross shows the location of a = 1.33 and 
b = —3.36, and the blue cross shows the location of a = 1.40 and 
b = —3.17, corresponding to the minimum of X for all values of 
a and b. 
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Age (Myr) Age (Myr) 

Fig. 5. Plots showing the influence of the free parameter a on the rotational evolution of solar mass stars starting at ~150 Myr 
at the 90th percentile of the rotational distribution. The left panel shows angular velocity against stellar age for different values of 
a, where the triangles show the observational constraints on the 90th percentile track given in Table [I] The horizontal dashed line 
shows the saturation threshold for both the mass loss rate and the magnetic field strength. The right panel shows corresponding 
wind mass loss rates along each track. 


classical Runge-Kutta method and Eqn. [2] to calculate the 
wind torque. 

The three free parameters that we fit later in this sec¬ 
tion are K T , a, and b. The parameter c in Eqn. [6] is impor¬ 
tant since it controls the mass dependence of the saturation 
threshold. We first fit K T and a by considering the example 
of the solar wind and the spin down of solar mass stars. We 
then fit c by first setting c = 2; this is close to the value sug¬ 
gested from the saturation of X-ray emission but does not 
allow us to get acceptable fits to the spin down of 0.5 M 0 
and 0.75 M 0 stars. Using this value, we find the best fit 
value for b. We then iteratively increase c, recalculating the 
best fit value of b each time until we get a solution that we 
judge by eye as being an acceptable fit. This allows us to get 
good fits to the rotational evolution models while keeping 
c as close to our expectations from X-ray observations. 

4.1. The current solar wind and the value of K T 

It is possible to constrain K T by considering the torque ex¬ 
erted on the Sun by the current solar wind. Observationally, 
the torque on the Sun averaged over long time periods can 
be easily estimated from observed spin down laws, such as 
Eqn. [?] The angular velocity of solar mass stars beyond 
~1 Gyr is known to evolve with age, t, as fl* oc t~ 0 ' 566 . 
This implies that the rate at which the angular velocity 
of the Sun is changing is dfl^/dt = — O.SGGf^tg 1 , where 
t 0 = 4.6 Gyr. Inserting this into Eqn.[IJ we get the torque 
from the current solar wind as 

r w>0 = — 0.5 (8) 

where J 0 = 6.87 x 10 53 g cm 2 . The first term on the RHS 
of this equation, involving d/ 0 /dt, is negligible. We get a 
torque exerted by the current solar wind on the Sun of 
—7.15 x 10 3 ° erg s _1 . This is the torque from the current so¬ 
lar wind averaged over time periods of hundreds of millions 
of years, and does not necessarily represent the solar wind 


torque at any instant in time. However, to calculate K T , 
we make the assumption that the current solar wind and 
magnetic field values are typical for the solar wind averaged 
over such long time periods. Inserting B^ip = 1.35 G and 
M* = 1.4 x 10 -14 M 0 yr -1 into Eqn.JSJgives a wind torque 
of —6.46 x 10 3 ° erg s _1 . Therefore, we find that K r = 11 
is necessary for our predicted wind torque to match the 
current solar wind. 


The fact that we need such a large value of K T could 
be interpreted in several ways. This could be a result of the 
fact that we only consider the dipole component of the field 
and ignore the other field components, which could cause 
us to underestimate the wind torque, especially at cycle 
maximum when the dipole field is negligible. Alternatively, 
it could be that we underestimate the values of Rdip,© and 


Mq . We derived these values using solar magnetic field and 
wind measurements that extend no more than a few decades 
into the past, whereas our determination of the solar wind 
torque using Eqn. [8] is sensitive to the values averaged over 
hundred of millions of years. If the Sun is currently at a long 
term level of low activity, then we would require a large 
value of K t to compensate for the low values of i?dip,o 
and Mq. Thi s interpretation, whil e highly speculative, is 
supported by Wright et al. (2011) who showed that the 
average solar X-ray luminosity is a factor of 2-3 below what 
would be predicted given the Sun’s mass and rotation rate. 


4.2. Fitting a and b 

In order to find the best fit values of a and b, we define a 
goodness-of-fit parameter that we call X. To calculate this 
parameter for a given set of values for the free parameters, 
we first calculate rotation tracks starting at the 10th, 50th, 
and 90th percentiles of the 150 Myr rotational distribution 
in each mass bin considered. For each of the observationally 
constrained percentiles listed in Table [T] we then calculate 
the square of the difference in the logarithm of the angular 
velocity between the observed value, H 0 b s , and the value 
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predicted by the model, f2 mo dei- We also include the pre¬ 
dicted rotation rates at 5 Gyr from Eqn. [7] for each mass 
bin, but only compare them to the 10 th percentile tracks. 
We do not attempt to encourage the 50th and 90th per¬ 
centile tracks to spin down to these values. Our goodness- 
of-fit parameter is then the sum of all of these values and 
is given by 


This is not exactly true since we have have ignored a 12* 
term in the denominator of Eqn. [2] but it is a good approxi¬ 
mation at these age^] Inserting M* oc f2“ and B^ip oc 12*- 32 , 
this becomes 


d!2* 

dt 


oc 12 °- 56 “+ 2 - 15 , 


( 12 ) 


X — ^ ( 7?. (log f2 0 bs,i lt)g 12 mo del,z) ? (9) 

i 

where the sum is the sum over all percentiles for each age 
and each mass bin, plus also the 5 Gyr predicted rotation 
rates. In this equation, 7 ,; is a weighting parameter that we 
use to increase or decrease the importance of certain parts 
of the observational constraints. For most of the observed 
percentiles from young clusters, we set 7 j = 1 , and for the 
comparison between the predictions of Eqn. for rotation 
at 5 Gyr and the 10th percentile models, we set % = 5. 
This is done to force the models to spin down to the desired 
rotation rates at later ages. Since one of the most difficult 
features of the observations to reproduce in the models is 
the slow spin down of rapidly rotating stars, especially at 
low masses, we set 7 * = 2 for the 90th percentiles of the 
distributions at 550 Myr. 

Other than K T , which we have already determined by 
considering the current solar wind, the only free parame¬ 
ter that influences the rotational evolution of solar mass 
stars is a. We therefore determine a by finding the value 
that minimises our goodness-of-fit parameter, X, consider¬ 
ing only the 1.0 M 0 mass bin. The value of a determines 
the dependence of the mass loss rate in the wind on stel¬ 
lar rotation. In the upper panel of Fig. [4] we show our 
goodness-of-fit parameter, X, against a and estimate that 
a ss 1.33. In Fig. [5] we show how the rotational evolution of 
the 90th percentile track, and the corresponding predicted 
mass loss rates, depends on our chosen value of a. Clearly 
both the predicted rotational evolution and the predicted 
mass loss rates are highly sensitive to the value of a. For 
low values of a, relative to the best fit value, a star starting 
at the 90th percentile at 150 Myr spins down too slowly 
and is well above the 90th percentiles at later ages. Simi¬ 
larly, high values of a lead to the star spinning down too 
fast and reaching the saturation threshold too early. Inter¬ 
estingly, the rotation rate of the different tracks at 5 Gyr is 
not highly sensitive to the value of a, though there is some 
influence. 

It is interesting to see if we can reproduce this value of a 
using different reasoning. At a given stellar mass, ignoring 
the change in stellar radius as the star ages, the mass loss 
rate is given by M* oc 12“. Since we know that at ages 
beyond 1 Gyr, the angular velocity has a dependence on 
age of 12 * oc t~ a - 5ee , we can estimate that 


d!2* 

dt 


oc 12 


2.77 

★ 


( 10 ) 


If we also ignore the evolution of the moment of inertia of 
the star, then dl 2 */dt is proportional to the wind torque, 
and therefore Eqn. [2] implies 


„ d0.87 7.V0.56 o 

-jf a B dip M * “*• 
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( 11 ) 


which compared with Eqn. |10| implies that a~ 1.1. This 
is similar to our determination of a ~ 1.33. Although this 
simple reasoning gives us a similar value of a to the more 
sophisticated analysis, it is still necessary to determine a 
by minimising our goodness-of-fit parameter because we 
are partly fitting a to the rotational evolution of rapidly 
rotating stars in the first billion years when the stars do 
not follow a Skumanich style spin down law. The fact that 
our determinations of a are so similar shows that the early 
rotational evolution of rapidly rotating solar-mass stars is 
perfectly consistent with the Skumanich style spin down at 
later ages. 

Now that we have derived the values of K T and a, we 
can find the best fit value of b. To do this, we consider the 
rotational evolution of all three mass bins. In the middle 
panel of Fig. [4j we show the dependence of X on b , assum¬ 
ing that a = 1.33. Our best fit model is the model where 
b « —3.36. 

Once again, it is interesting to see if we can estimate a 
similar value of b using different reasoning. From Eqn. [2j 
the wind torque on a star is approximately given by 


r w oc £°:fM°' 56 i? 2 - 87 ft*. ( 13 ) 

When determining K T , we showed that dtt+fdt oc 
which naturally follows from assuming a power-law depen¬ 
dence of 12* on t. Ignoring the time evolution of the moment 
of inertia, this implies that at a given age, r w oc 1*12*. Since 
the mass loss rate is given by M* oc this becomes 

1*12* oc B^ 7 Rl "nl- 74 M°- 56b , (14) 

where we have inserted a = 1.33. We can see from this 
analysis that our best fit value of b is likely to depend, at 
least weakly, on our best fit value of a. We can estimate 
the value of b by considering the ratio of the torques on 
two stars that have the same age. If we consider two stars, 
Eqn. [T4| implies that 



(15) 

4 Ignoring the denominator in Eqn. [ 2 ] is reasonable for 
slow rotation. For example, the current solar parameters give 
K%Ve sc = 9.8 x 10 12 cm 2 s -2 and = 3.5 x 10 10 cm 2 s -2 , 

meaning that the rotational dependence is negligible. At fast 
rotation, the 12*_R* term dominates the A' 2 u 2 ac term and so the 
approximation given in Eqn. m is not reasonable. For 0.5 Mg 
and 1.0 Mq, the two terms become equal at rotation rates of ap¬ 
proximately 30f2Q and 1712g respectively. We take into account 
all terms in Eqn. [2] in our derivation of a « 1.33 and b « —3.36. 
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where quantities marked with a prime are for one star and 
quantities not marked with a prime are for the other star. 
Rearranging for b gives 


log 


b = 


©(&)'“’ (£)"“(£) 


-3.99 


-0.74' 


0.56 


lo S (£) 


(16) 


Consider two 5 Gyr old stars with masses of 0.5 M 0 and 
1.0 Mq. From Eqn. 5] and Eqn. [7 we can estimate the 
equatorial dipole fielclstrengths ana rotation periods, and 
from the stellar evolution models discussed in Section [2] we 
can estimate the stellar radii and moments of inertia. For 
the 0.5 Mq star, we estimate a rotation period of ~51 days, 
a dipole field strength of 1.85 G, a radius of 0.53 Rq, and a 
moment of inertia of 2.1 x 10 53 g cm 2 . For the 1.0 M 0 star, 
we estimate a rotation period of ~27 days, a dipole field 
strength of 1.25 G, a radius of 1.03 R 0 , and a moment of 
inertia of 6.9 x 10 53 g cm 2 . Inserting these quantities into 
Eqn. [l6]gives b ~ —4.1. This is similar to our derived value 
of b « —3.36. 

To check that our estimates of a and b are robust, we 
fit the values once more by minimising the goodness-of-fit 
parameter, X. Above, we fit the two parameters separately, 
but to check that this gives us approximately the best fit, 
we fit the two parameters simultaneously by producing a 
grid of rotational evolution models with ranges of values of 
a and b. In the lower panel of Fig. [4] we show a contour plot 
of X from this grid. Clearly there is some degeneracy be¬ 
tween our values of a and b , which can be understood from 
our above analysis. The value of a is well constrained, with 
1.2 and 1.6 likely representing reasonable limits. The con¬ 
straints on b are much looser, with reasonable lower and 
upper limits being approximately -4 and -2 respectively. 
The red cross shows our best fit of a ss 1.33 and b ss —3.36, 
which is clearly very close to the minimum of A", which lies 
at a « 1.40 and b ss —3.17. We choose to use the first set 
of values of a and b since they give the best fit to the ro¬ 
tational evolution of solar mass stars, which is much better 
constrained beyond 1 Gyr observationally than the rota¬ 
tional evolution of lower mass stars. 

One interesting result of these fits is the fact that at 
a given rotation rate in the unsaturated regime, low-mass 
stars have higher mass loss rates per unit surface area than 
high-mass stars. We predict that a 0.5 M 0 star has a mass 
loss rate per unit surface area that is a factor of ten larger 
than a 1.0 M 0 star with the same rotation rate. This trans¬ 
lates into a mass loss rate that is a factor of ~3.4 higher for 
the 0.5 M 0 star. That low-mass stars have higher mass loss 
rates is an interesting result. Assuming that R* oc M 3 8 , 
this is inevitable if b < —1.6, which Fig. [4] shows clearly 
must be the case for our rotational evolution models to 
provide acceptable fits to the observational constraint^] 


5 Recently, Gallet & Bouvier (2015) produced rotation mod¬ 
els for stars with masses of 0.5 Mq and 0.8 M 0 . To prescribe 
the ma ss loss rates, they used the model of |Cranmer fe Saar| 
(201 lb|), which predicts much lower values for low mass stars 
(for the 0.5 M 0 mass bin, the mass loss rate is between 10 1 ' 
and 10- 16 M 0 yr -1 ), in strong disagreement with our results. 
However, to re produce the observational constraints, |Gallet fc| 
Bouvier (20151 had to artificially increase the multiplicative con¬ 
stant in their formula for the wind torque for low mass stars, 


Vidotto et al. (2014b) presented 3D MHD models of 
the winds from six Mndwarfs based on realistic non-dipolar 
magnetic field geometries measured using the Zeeman- 
Dopple r Imaging technique. In Table 3 of |Vidotto et al.] 
(2014bI, power-laws for the relations between different 
quantities were given based on fits to the results of their 
models. They found that the wind torque is related to the 
mass loss rate as r w oc M 2A8± °- 56 . This might appear to be 
in contradiction to the M 356 dependence derived by Matt 
et al. (2012| ) that we use in this paper, but it is important to 
recognise tnat the relation of Vidotto et al. (2014b I is not 
the dependence of wind torque on mass loss rate, but is the 
correlation between wind torque and mass loss rate, and 
therefore also includes the correlations between mass loss 
rate and other quantities, such as magnetic field strength 
and rotation rate. At a g iven stellar m ass and radius, the 
torque formula derived by Matt et al. (2012) approximately 
reduces to 


oc R° ; 87 Ae 56 fi* 


Our derivation of the de¬ 
pendence of mass loss rate on rotation and Eqn. 0 imply 
that fl* oc M^ 1 ' 33 and i?dip oc AV* ' 32 ^ 133 . Therefore, we 
find th at r w oc M 218 , in exce llent agreement with the re¬ 
sults of Vidotto et al. (2014b] (the apparent perfect agree¬ 
ment is of course coincidental). This agreement suggests 
that not considering non-dipolar magnetic field geometries 
in our rotational evolution model does not influence our 
results significantly. 

It is possible with our results to construct simple scaling 
laws for the wind torque, r w , as a function of stellar rotation 
and mass. Assuming R* oc Af 9 - 8 means that the mass loss 
rates and dipole field strengths vary as M* a fi* 33 M “ 1,76 
and Bdip oc fi* 32 M* * 43 , where the latter equation was de¬ 


rived by assuming a M~ lm dependence to the convective 
turnover time dReiners et al.ll20T4l) . These are true until 
the saturation rotation rate given by fl sat oc M 2 ' 3 . Insert¬ 
ing these scaling laws into the simplified version of the wind 
torque formula given in Eqn.[l3]gives r w in the unsaturated 
regime as 


T w 


^w,© 



(17) 


where r w 0 ss —7.15 x 10 3 ° erg s 1 is the current solar 
wind torque. In the saturated regime, we find instead that 


7”w 


T w,© 


15 1 ' 89 



\MrJ 


(18) 


Given that Eqn. [13] is only an accurate approximation at 
slow rotation, the rotation dependence in this formula is 
not reasonable for the fastest rotators. Interestingly, the 
wind torque has different dependences on stellar mass in 
the saturated and unsaturated regimes. This is due simply 
to the strong mass dependence of the saturation threshold. 


which has a similar effect as the M b term in our Eqn 4. Our re- 
su lts correspond to the int erpretatio n, discussed in Section 5.3 
of Gallet & Bouvier (20151, that the Cranmer & Saar (201 lb I 
model underestimates the mass loss rates of low mass stars. IGal-l 


let & Bouvier (20151 suggested that the mass loss rates of the 
0.5 Mq stars should be increased by a factor of 860 relative to 
what they would otherwise predict, which would be much more 
consistent with our results. 
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This explains the interesting fact that rapidly rotating low- 
mass stars spin down slower than rapidly rotating high- 
mass stars, but slowly rotating low-mass stars spin down 
quicker than slowly rotating high-mass stars. 

4.3. Sources of uncertainty in the fit parameters 

We now have best fit values for the free parameters in our 
rotational evolution model. Constraining these free param¬ 
eters is important because it not only allows us to predict 
the rotational evolution of stars on the main-sequence, but 
gives us an observationally constrained scaling law for mass 
loss rate as a function of stellar mass, radius, and rota¬ 
tion. In the stellar wind model developed in Paper I, the 
largest unknown was how to scale the base temperature 
and density to other stars. This scaling law for the mass 
loss rate gives us a large part of the solution to this prob¬ 
lem. However, before we analyse what these results mean 
for the evolution of wind properties on the main-sequence, 
we should discuss sources of uncertainty in our results. We 
concentrate mostly on the parameters a and b since they 
are the most important for our wind model. Probably the 
main source of uncertainty in our predicted values of a and 
b comes from how we relate stellar and wind properties to 
wind torques. Although our results are also based on stel¬ 
lar evolution models and observed rotation rates of stars at 
different ages, these are much better constrained. We em¬ 
phasise that we are only speculating in this section about 
reasons why our results might be inaccurate and have no 
clear reason to think that there are such inaccuracies. 

The driver of our rotational evoluti on model is the 
torque formula (Eqn. [2]) derived by Matt et al. (20121, which 


Vidotto et al. 

2009) 

perature in al 

mod 


Matt et al. 


(2012) also set the tem- 


of Matt et al. 


would influence the models in a similar way as a and b, i.e. 
both assumptions add power law dependences on and 
M* into the wind torque. 

Another source of uncertainty is in how we calculate the 
strength of the dipole component of the field for each star 
given its mass and rotation rate. We do this using Eqn. [5j 
which is based on the correlation between the magnetic field 
strength and rotation derived by Vidotto et al. (2014a). The 
form of this relation is B^ip oc Ro d . The exact value of d, 
howeve r, is uncertain since the relation of |Vidotto et al.| 
(2014a I is based on measurements of the large scale com¬ 


allows us to calculate wind torques from the stellar mass, 
radius, rotation rate, dipole magnetic field strength, and 
mass loss rate. This equation is derived from 2D MHD wind 


models. Matt et al. (2012) assumed a polytropic equation 
of state, with a = 1.05 everywhere, and a base temperature 
and density of the wind that is uniform over the stellar sur¬ 
face. Such models produce realistic distributions of wind 
speeds for a given magnetic field geometry, but are inaccu¬ 
rate in that the contrast between the slow and fast wind 
is underestimated (for example, s ee Fig . 5 and Fig. 6 of 


is a fixed fraction of the escape velocity at the base of the 
wind. Finally, they assumed that the stellar magnetic fields 
are all dipolar. In real winds, these assumptions might not 
be fulfilled, which could lead to the dependences of torque 
on the various parameters being different to the predictions 

( |20l2 l. 


We have needed to add an extra free parameter, K T , 
into our model for calculating the torque (i.e. Eqn. [3]). The 
free parameter K T is not a meaningless fudge factor, but 
represents real physical processes that we do not under¬ 
stand, such as the details of the wind driving mechanisms 
and non-dipolar field geometries. Also, I\ T is influenced by 
our measurements of the current solar mass loss rate and 
dipole field strength, which might not be fully represen¬ 
tative of the values averaged over long time periods. We 
assume that K T is the same for all stars; however, it could 
be that this parameter should also have a dependence on 
stellar mass and rotation. We could take this into account 
by assuming K T oc Cl^M{ and using e and / as free param¬ 
eters in the rotational evolution models. However, e and / 


ponents of the magnetic field strength, which is not exactly 
the same as the dipole component of the field, and since the 
relation between the magnetic field strength and the Rossby 
number contains a lot of scatter, from which a range of val¬ 
ues of d can provide acceptable fits. The parameter d has the 
same effect on the rotational evolution models as a, since 
they both influence the power law dependence of the wind 
torque on angular velocity. A change in d would lead to a 
change in the opposite direction in the best fit value of a. 
In addition, it is not completely clear that the dipole com¬ 
ponent of the field should be scaled with Rossby number 
instead of some other parameter. Using the Rossby number 
adds a mass dependence into the dipole component of the 
field, such that at a given rotation rate, lower mass stars 
have larger magnetic field strengths than higher mass stars. 
This dependence has a similar effect as the M * dependence 
in our assumption about the mass loss rate. If instead, we 
scaled Rdip with f2*, we would then need a stronger mass 
dependence for the mass loss rate to compensate in the ro¬ 
tational evolution models. 

Looking at the above discussion in a more general way, 
by making the assumption that Af* oc we are 

adding extra power law dependences on rotation and mass 
into the wind torque which we then interpret as being en¬ 
tirely due to variations in the mass loss rate. If there were 
other sources of such dependences then we would confuse 
them with the dependence of mass loss rate on stellar ro¬ 
tation and mass. 

In the previous section, we are able to estimate values of 
a and b that are similar to our best fit values using cruder 
analytic reasoning. This is both encouraging and discour¬ 
aging. It is encouraging because one of the major sources 
of uncertainty in our model is the angular velocity at which 
saturation occurs for the mass loss rate and dipole field 
strength. Since the cruder reasoning only considers rota¬ 
tion at later ages, we can see that this uncertainty does not 
significantly influence our result. It is discouraging because 
it means that our best fit value of b is mostly sensitive to 
the spin down of low mass stars at later ages, which is con¬ 
strained by Eqn. [7] However, the spin down of stars at later 
ages, as predicted by Eqn.[^ is constrained almost entirely 
by the Sun, and so additional observational confirmation of 
the spin down of low-mass stars is desirable. 

Another source of uncertainty in our model comes from 
the possibility that young stars arrive at the ZAMs with 
inner cores that are rotating faster than the stellar surfaces. 
Core-envelope decoupling would lead to a spin up torque 
acting on the stellar surface that would work against the 
wind torque, meaning that in order to reproduce the same 
rotation track, a larger spin down torque, and therefore a 
larger mass loss rate, would be needed from the wind. For 
solar mass stars, this w ould lead to our predicte d value of 
a being larger. However, Gallet & Bouvier (20131 predicted 
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from rotational evolution models of solar mass stars that 
around the ZAMS, the core would be rotating a factor of 
~2 faster than the envelope, and that by the age of the 
Sun, the core and the envelope would have the same angular 
velocities. In this scenario, our wind would need to remove 
twice as much angular momentum from the star between 
the ZAMS and 5 Gyr, which is unlikely to require a much 
larger value of a. 

In our models, we consider only the quiet winds of stars. 
It has recently been suggested that the winds of more active 
stars could be dominated by coronal mass ejections (CMEs; 
Khodachenko et al.||2007| |Aarnio et al.|[20l2| |Drake et al.| 


2013). Lacking more detailed physical knowledge of CME 


ejections from more active stars, the assumption that the 
angular momentum removal per unit mass loss is the same 
for hypothetical CMEs as for quiet winds implies that our 
scaling law for the mass loss rate remains valid. Should a 
CME remove more (less) angular momentum per unit mass 
loss than the quiet wind, we would overestimate (under¬ 
estimate) the wind mass loss rate when fitting our scaling 
law for M to the observations (if CMEs do indeed carry 
away a significant amount of mass and angular momen¬ 
tum). Given the fundamental uncertainties of these ques¬ 
tions and the lack of observation al support for strong C ME 
dominated winds in active stars (Leitzinger et al. |2014[ ), we 
do not consider the question of CME angular momentum 
transport further. 


5. Rotational Evolution: Results 


In Fig. § we show the rotational evolution tracks for our 
best fit models. The three panels show the evolution of stars 
with masses of 0.50 M 0 , 0.75 M 0 , and 1.00 M 0 . The obser¬ 
vational constraints are shown as triangles for the 10th and 
90th percentiles and squares for the 50th percentiles which 
are summarised in Section [3] and Table [T] The squares at 
older ages are predictions of the 5 Gyr rotation rates from 
the gyrochronol ogical relation derived by |Mamajek fe Hil-| 
lenbrand (20081. Our model clearly fits the observational 


constraints very well, especially in the 0.5 M 0 and 1.0 M 0 
mass bins. 

For solar mass stars, we are able to explain the rota¬ 
tional evolution of stars from 100 Myr to 5 Gyr well. This 
is done with only two free parameters, K T and a. According 
to our model, a star at the 90th percentile of the rotational 
distribution at 150 Myr will remain in the saturated regime 
until about 300 Myr. By 500 Myr, the tracks have mostly 
converged, and by 1 Gyr, there is very little difference be¬ 
tween the tracks. Our tracks probably overestimate slightly 
the spread in rotation rates at 1 Gyr compared to the ob¬ 
servational constraints, though the difference is very small. 

It is interesting t o compare our models to those of |Gal-| 
let & Bouvier (20131 for solar mass stars. The main differ- 
ence between our model and the model of Gallet & Bo uvierl 


(20131 on the main-sequence is the use of core-envelope de¬ 
coupling. Gallet & Bouvier (2013) found that without the 
extra spin-up torque from the core on the envelope due to 
core-envelope decoupling, rapidly rotating stars would spin 
down too quickly by the age of ~650 Myr. As we show in the 
next section, our model for M t provi des similar results to 
the model of Cranmer & Saar (2011b I used in their model, 


so the difference is unlikely due to difference values of M+. 
The difference is likely that we use a weaker scaling for B^ip 
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Fig. 6. Plots showing the rotational evolution of stars at the 
10th, 50th, and 90th percentiles of the rotational distribution 
from 100 Myr to 5 Gyr. The three panels show the rotation 
tracks for stars with masses of 0.50 Mg ( upper panel), 0.75 M 0 
(middle panel), and 1.00 M 0 ( lower panel). The observational 
constraints on the rotational evolution of stars at the 10th per¬ 
centile (downward pointing triangles), 50th percentile (squares), 
and 90th percentile (upward pointing triangles) of the rotation 
distribution. The horizontal dashed lines show the saturation 
thresholds for the mass loss rates and dipole field strengths. 
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Fig. 7. Figure showing how rotation in our composite cluster evolves between 100 Myr and 2 Gyr according to our rotational 
evolution model. The cluster is made up of ~1500 stars in the four young clusters Pleiades, M50, M35, and NGC 2516, all of which 
have ages ~150 Myr. The upper row shows the distribution predicted by evolving this composite cluster back to 100 Myr, and the 
other three rows show the distribution evolved forward to 500 Myr, 1 Gyr, and 2 Gyr. The left column shows the position of each 
star in the distribution and the right column shows the 10th, 50th, and 90th percentiles of the distribution as a function of mass 
for each age. The dashed lines show the saturation threshold determined by Eqn. [6] 
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Stellar Mass (M 0 ) 

(fi 0 ) 

(0©) 

^10 (^©) 

^50 (^©) 

^90 (^©) 

Mio (M e ) 

-^50 (M©) 

Mgo (Mq) 

100 Myr 

All 

27.43 

44.62 

3.24 

8.49 

78.79 

- 

- 

- 

0.5 

32.39 

35.09 

3.88 

18.05 

85.25 

11.84 

11.84 

11.84 

0.6 

29.04 

35.28 

3.12 

11.13 

88.65 

8.49 

14.36 

14.36 

0.7 

27.22 

42.93 

3.02 

6.52 

82.76 

6.26 

17.42 

17.72 

0.8 

25.44 

52.11 

3.07 

5.20 

69.82 

5.08 

10.21 

21.13 

0.9 

22.65 

54.66 

3.23 

5.64 

53.11 

4.50 

9.48 

25.20 

1.0 

23.71 

51.17 

3.76 

7.15 

60.10 

4.45 

10.46 

28.02 

250 Myr 

All 

16.43 

28.08 

2.67 

5.19 

50.29 

- 

- 

- 

0.5 

23.00 

26.70 

2.73 

12.11 

62.18 

10.29 

11.88 

11.88 

0.6 

19.07 

24.65 

2.51 

6.35 

60.91 

6.37 

14.38 

14.38 

0.7 

15.62 

27.11 

2.56 

3.89 

45.81 

5.02 

8.78 

17.76 

0.8 

13.56 

32.32 

2.62 

3.77 

34.54 

4.11 

6.67 

21.18 

0.9 

11.07 

30.92 

2.82 

4.24 

16.39 

3.78 

6.52 

25.30 

1.0 

10.74 

23.39 

3.25 

5.03 

14.74 

3.74 

6.69 

27.94 

500 Myr 

All 

7.56 

12.68 

2.11 

3.28 

19.24 

- 

- 

- 

0.5 

12.47 

15.90 

1.96 

5.77 

33.83 

6.63 

11.94 

11.94 

0.6 

9.11 

12.26 

1.96 

3.04 

28.13 

4.58 

8.21 

14.41 

0.7 

6.05 

11.03 

2.08 

2.69 

13.93 

3.84 

5.39 

17.82 

0.8 

5.23 

12.41 

2.19 

2.82 

6.82 

3.25 

4.57 

14.75 

0.9 

4.47 

9.67 

2.37 

3.16 

4.77 

3.02 

4.43 

7.66 

1.0 

4.12 

2.88 

2.69 

3.64 

4.98 

3.00 

4.49 

6.81 

750 Myr 

All 

4.10 

5.82 

1.74 

2.52 

7.43 

- 

- 

- 

0.5 

6.60 

8.87 

1.57 

2.84 

17.12 

4.99 

10.96 

12.00 

0.6 

4.48 

5.45 

1.61 

2.16 

11.89 

3.55 

5.22 

14.43 

0.7 

3.00 

3.73 

1.79 

2.17 

4.34 

3.15 

4.07 

10.22 

0.8 

2.84 

3.56 

1.91 

2.37 

3.55 

2.73 

3.62 

6.21 

0.9 

2.82 

2.00 

2.08 

2.62 

3.34 

2.55 

3.48 

4.80 

1.0 

2.98 

0.61 

2.36 

2.93 

3.57 

2.59 

3.45 

4.50 

1000 Myr 

All 

2.66 

2.82 

1.48 

2.08 

3.49 

- 

- 

- 

0.5 

3.60 

4.74 

1.32 

1.89 

8.33 

3.96 

6.41 

12.06 

0.6 

2.57 

2.19 

1.41 

1.78 

4.70 

2.97 

4.05 

14.46 

0.7 

2.11 

0.99 

1.58 

1.88 

2.82 

2.67 

3.36 

5.77 

0.8 

2.16 

0.81 

1.71 

2.06 

2.67 

2.36 

3.01 

4.28 

0.9 

2.29 

0.45 

1.88 

2.27 

2.71 

2.25 

2.88 

3.66 

1.0 

2.51 

0.37 

2.10 

2.50 

2.93 

2.28 

2.86 

3.55 

1500 Myr 

All 

1.66 

0.72 

1.16 

1.60 

2.08 

- 

- 

- 

0.5 

1.56 

1.21 

1.00 

1.27 

2.22 

2.75 

3.79 

7.98 

0.6 

1.48 

0.36 

1.15 

1.38 

1.96 

2.28 

2.90 

4.63 

0.7 

1.56 

0.26 

1.31 

1.53 

1.91 

2.09 

2.57 

3.45 

0.8 

1.68 

0.24 

1.43 

1.67 

1.95 

1.87 

2.30 

2.83 

0.9 

1.82 

0.22 

1.59 

1.83 

2.07 

1.82 

2.20 

2.58 

1.0 

1.98 

0.20 

1.76 

1.98 

2.22 

1.88 

2.19 

2.56 

2000 Myr 

All 

1.33 

0.30 

0.96 

1.34 

1.67 

- 

- 

- 

0.5 

1.09 

0.32 

0.82 

1.04 

1.41 

2.15 

2.93 

4.39 

0.6 

1.19 

0.18 

0.98 

1.16 

1.45 

1.85 

2.32 

3.10 

0.7 

1.32 

0.17 

1.13 

1.31 

1.53 

1.73 

2.10 

2.60 

0.8 

1.43 

0.16 

1.25 

1.44 

1.61 

1.58 

1.90 

2.21 

0.9 

1.56 

0.15 

1.39 

1.57 

1.72 

1.54 

1.81 

2.04 

1.0 

1.67 

0.13 

1.54 

1.67 

1.82 

1.62 

1.81 

2.03 


Table 2. Table giving statistical properties of our composite cluster at different ages. The cluster is composed of ~1500 stars with 
known rotation periods and masses from Pleiades, M50, M35, and NGC 2516, all of which have ages ~150 Myr. We have evolved 
the cluster backwards to 100 Myr and forwards to 2 Gyr using our rotational evolution model. The distributions of rotation rates 
at ages of 100 Myr, 500 Myr, 1 Gyr, and 2 Gyr for this cluster are shown in Fig. [7] In this table, we give, from left to right, the 
average angular velocity, fi*, the standard deviation in the angular velocity, <rn, the 10th, 50th, and 90th percentiles in the angular 
velocity distributions, and our predicted wind mass loss rates at these percentiles. For each age, we give the values for all of the 
stars between 0.4 Mg and 1.1 Mg, and for masses of 0.4 Mg, 0.5 Mg, 0.6 Mg, 0.7 Mg, 0.8 Mg, 0.9 Mg, and 1.0 Mg. At each 
specific mass, we calculate the values of each parameter using all stars within 0.1 Mg of the given mass. 
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with rotation rate. In our mode l -Bdip oc fl 3 ' 32 , whereas in 
the model of Gallet & Bouvier (20131, a stronger scaling 
closer to Bdip oc was used. This means that we pre¬ 
dict weaker torques for rapidly rotating stars, allowing them 
to spin down slower without core-envelope decoupling. We 
warn however that we have only shown that the evolution 
of solar mass stars on the main-sequence can be explained 
without the use of core-envelope decoupling. However, core¬ 
envelope decoupling is also an important factor on the PMS, 
which is a phase of rotational evolution that we do not con¬ 
sider. 

The rotation tracks for the lowest mass stars that we 
consider are shown in the upper panel of Fig. [6] Our models 
clearly provide a good fit to the observational constraints, 
though there is some indication that we slightly underes¬ 
timate the spread in rotation rates at 550 Myr. The dif¬ 
ference, however, is very small. There are two main differ¬ 
ences between the rotational evolution of 1.0 M 0 stars and 
0.5 M 0 stars. Firstly, in the first Gyr of evolution just af¬ 
ter the ZAMS, much less spin down of the rapidly rotating 
stars takes place. Secondly, low-mass slowly rotating stars 
spin down faster than higher mass stars. This results in the 
slowly rotating tracks in young clusters having a shape such 
that higher mass stars rotate faster than lower mass stars, 
with a similar trend at later ages. 

The most noticeable difference between the models and 
the observations are in the 0.75 M 0 mass bin. All three 
tracks spin down to approximately the same rotation rate 
by 5 Gyr, but this value is slightly larger than the pre¬ 
dicted value. It is unclear if this is a result of the model 
underestimating the wind torques at later ages for 0.75 M 0 
stars, or if Eqn. [7] overestimates the rotation periods for 
stars of this mass. This latter suggestion is plausible given 
that the exact form of the rotation period dependence on 
stellar mass at a given age beyond 1 Gyr is not well con¬ 
strained observationally, but is the result of assuming that 
the shape of the slowly rotating tracks seen in young clus¬ 
ters remains the same as the stars evolve to later ages. The 
nearby ~0.7 M 0 star 61 Cyg A provides a test for which of 
the tw o predictions is likely to be more realistic. [Donahue] 


(19961 esti mated that 61 Cyg A has a rotation period 
So days and Kervella et al. (2008) estimated an age for 


et al. 
of 

the star of 6.0 ±1.0 Gyr. This suggests that our prediction 
for the s pin down of stars with this mass are more likely. 
However, Donahue et al. (1996) also estimated a similar ro¬ 
tation period of ~38 days for the lower mass companion, 
61 Cyg B, which is difficult to reconcile with both our pre- 
dictions and the g yrochronological relation of |Mamajek fc| 
Hillenbrand (20081 if the age of 6 Gyr is reasonable. 

In Section [3] we combine the rotation period measure¬ 
ments for the Pleiades, M50, M35, and NGC 2516. These 
are four young clusters with ages of ~150 Myr, and we 
assume that the combination of the rotation periods for 
the four clusters represents the rotational distribution at 
150 Myr. To illustrate the results of our rotational evolution 
model, we evolve this cluster back to 100 Myr, and forward 
to 2 Gyr. In Fig. [7] we show the distribution at 100 Myr, 
500 Myr, 1 Gyr, and 2 Gyr. The left column shows the 
locations of each star in the distribution at each age and 
the right column shows the 10th, 50th, and 90th percentiles 
of the distributions as a function of mass. The percentiles 
were calculated at each mass by considering all stars within 
0.1 M 0 of that mass. In Table [2] we give these percentiles 
for masses of 0.5 M 0 , 0.6 M 0 , 0.7 M 0 , 0.8 M 0 , 0.9 M 0 , and 
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Fig. 8. Plot showing the fraction of stars in the saturated 
regime as a function of age based on the results shown in Fig. [7] 
for different mass bins. The solid lines correspond to the val¬ 
ues calculated using our constraints on the mass dependence of 
the saturation threshold given by Eqn. [6] with c = 2.3, and the 
dotted lines show the values calculated using the empirical sat¬ 
uration threshold for X-ray emission derived by |Wright et al. | 
( 20111 . 


1.0 M 0 at ages of 100 Myr, 250 Myr, 500 Myr, 750 Myr, 
1 Gyr, 1.5 Gyr, and 2 Gyr. As before, we assume mass bins 
with widths of 0.2 M 0 . The dashed lines in Fig. [7]show the 
mass-dependent saturation threshold. 

Fig- 03 shows clearly how the rotation rates of stars in 
each mass bin converge as they age. At 100 Myr, the 10th 
and 90th percentiles are almost independent of mass, indi¬ 
cating that the slowest and fastest rotators at each mass 
have similar rotation rates. However, the 50th percentile is 
strongly mass dependent. At all masses, a significant frac¬ 
tion of stars are in the saturated regime, especially for low 
mass stars where almost all stars are saturated. By 500 Myr, 
the 90th percentile has converged for solar mass stars, but 
is still at very fast rotation for low-mass stars. Almost all 
stars of solar mass have evolved out of the saturated regime, 
though a significant fraction of stars at low masses remain 
saturated. A similar distribution is seen at 1 Gyr, though 
the 90th percentile track has made more progress in con¬ 
verging with the other tracks at low masses and very few 
stars at low masses are saturated. Between 1 Gyr and 2 Gyr, 
the remaining rapid rotators converge, and at 2 Gyr, all 
stars lie on one track, with lower mass stars rotating slower 
than higher mass stars. At 2 Gyr, there are no remaining 
saturated stars in the distribution. The age dependence of 
the fraction of stars in the saturated regime for the different 
mass bins is shown in Fig. [8] 



6. The Evolution of Stellar Winds on the 
Main-Sequence 

It is now possible to couple the stellar wind model devel¬ 
oped in Paper I with the rotational evolution model that we 
develop in this paper to predict how stellar wind properties 
evolve on the main-sequence for a range of stellar masses. 
We review our wind model in Section |6.1| and then apply 
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it to the evolution of the solar wind in Section 16.21 and to 
stars with a range of masses in Section [673] 


6.1. Stellar wind model 


In Paper I, we construct a stellar wind model based on scal¬ 
ing the solar wind to other stars. Our wind model was calcu¬ 
lated using a ID numerical hydrodynamic model of the solar 


wind run using the Versatile Advection Code (Toth 1996 


Toth 


1997). The fundamental driving mechanism of our 


wind is thermal pressure gradients. In order for a thermal 
pressure driven wind to gain enough energy to accelerate to 
the correct wind speeds, it is necessary to heat the wind as 
it expands. To do this, we assume that pressure and density 
are related by a polytropic equation of state, p oc p a . In or¬ 
der to accurately describe the wind, we choose a value of a 
that varies with radial distance from the solar surface. Once 
a is set, the two free parameters that determine the wind 
properties are the base temperature, To, which determines 
the wind acceleration and strongly influences the wind den¬ 
sities, and the base density, no, which strongly influences 
the wind densities only. 


In Paper I, we run a grid of 1200 numerical wind models 
with a range of stellar masses, stellar radii, and wind tem¬ 
peratures. We show that the wind speed can be predicted 
from Tq only with the equation 


uiau « 73.39 + 224.14Xo - 11.28T 0 2 + 0.28T o 3 , (19) 


where uiau is the wind speed at 1 AU in km s 1 and T 0 
is in MK. The dependence of wind speed on To is shown in 
Fig-0 Once To is known, the base density can be calculated 
in order to give the required densities far from the star. 


For the solar wind, we constrain To and no using space¬ 
craft measurements. The solar wind breaks down quite 
clearly into two components based on the speed of prop¬ 
agation. The spatial distribution of slow and fast wind 
is closely a ssociated with the struc t ure of the solar mag - 
netic field (Wang & Sheeley 1990 Arge & Pizzo 20001. 
The slow wind travels approximately 400 km s” 1 and at 
1 AU has typical proton densities and temperatures of 
~ 5 cm -3 and 7.5 x 10 4 K. The fast wind travels approxi¬ 
mately 760 km s” 1 and at 1 AU has typical proton densities 
and temperatures of ~ 2 cm” 3 and 1.8 x 10 5 K. Fitting our 
models to the spacecraft measurements, we found base tem¬ 
peratures of 1.8 MK and 3.8 MK for the slow and fast com¬ 
ponents of the solar wind respectively. Our models provide 
an excellent description of the real solar wind far from the 
solar surface, but are unrealistic within the solar corona. 
In our model, we assume that the winds of all low-mass 
main-sequence stars break down into slow and fast wind 
components, and that the two components scale to other 
stars in the same way. 


For stars other than the Sun, spacecraft measurements 
are not available. We use instead the constraints on the 
mass loss rates from our rotational evolution model to de¬ 
rive the wind mass fluxes, pv. However, disentangling p and 
v is difficult and requires that we know the wind tempera¬ 
ture. In Paper I, we develop two models for scaling the base 
temperature. 



radius (AU) 



Fig. 9. Figure demonstrating the dependence of wind speed on 
the wind base temperature in our model. The upper panel shows 
wind speed against radial distance from the stellar surface for a 
solar mass and radius star for five different base temperatures. 
The lower panel shows the dependence of wind speed at 1 AU 
on To, as given by Eqn. |19| 


In Model A, we follow Holzwarth & Jardine (2007) and 
assume that the base temperature is proportional to the 
coronal temperature, such that 


To — T 0 ,q 


T„ 


1 cor,© 


( 20 ) 


where T cor is the plasma temperature averaged over the 
entire corona. In order to do this, we use the result of 
Johnstone & Giidel (2015) that there exists one univer¬ 
sal relation between coronal average temperature, T cor , 
and X-ray surface flux, Fx = Tx/(47ri? 2 ), for all low- 
mass main-sequence stars, such that T cor ss O.llT^ 26 , 
where T cor is in MK and Fx is in erg s 


cor 

_1 cm” 2 . We use 


20111 to esti- 


the empirical relations of Wright et al. 
mate Fx as a function of stellar mass and rotation. 


In Model B, we follow Matt et al. (2012) and assume 
that the sound speed at the base of the wind is a fixed 
fraction of the escape velocity, such that 


To = 


2 Gpm p 
7^B 


c s 

^esc 


M* 

R* 


( 21 ) 
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where the fraction c s /v esc is a constant that we con¬ 
strain from the solar wind models For the slow and 
fast winds respectively, we find c s /v esc = 0.329 and 
c s /v esc = 0.478. 

The two models are by definition equivalent for the Sun, 
but for other stars they can lead to different results. 

The final effect that we consider in our wind model is the 
acceleration of the wind from magneto-rotational effects. As 
an ionised wind expands away from a rotating magnetised 
star, it gains kinetic energy and angular momentum from 
the magnetic field, which can lead to increased wind speeds 
far from the star and, in extreme cases, increased mass loss 
rates. Due to this acceleration, for rapidly rotating stars, 
the wind structure becomes highly latitude dependent and 
properly taking this into account would require 2D or 3D 
MHD models which would be inappropriate in this paper. 
On the other hand, the influence of magneto-rotational ef¬ 
fects on the wind speeds is likely to be so significant for the 
most rapidly rotating stars that it would also be inappro¬ 
priate to ignore them entirely. For simplicity, we assume the 
influence of the rotation of the star can be charactarised by 
the Michel velocity, given by 


um = 


M* 


( 22 ) 


As we justify in Paper I, £?* = 0.2Bdi P , where i?dip is the po¬ 
lar strength of the dipole component of the magnetic field. 
In the following section, we discuss our model for calculat¬ 
ing Bdip for other stars. We follow Belcher & MacGregor 
([1976) and consider two regimes: these are the slow mag- 
netic rotator (SMR) and the fast magnetic rotator (FMR) 
regimes. A star is in the SMR regime when the speed that 
the wind would have in the absence of the star’s rotation 
is much larger than dm, in which case, magneto-rotational 
effects can be ignored. A star is in the FMR regime when 
wm is larger than the speed that the wind would have in the 
absence of rotation, in which case, the terminal velocity of 
the wind is approximately um ■ The current solar wind has a 
Michel velocity of ~40 km s -1 , and is therefore in the SMR 
regime. While this model is crude, it gives us a simple way 
to estimate the influence of magneto-rotational acceleration 
of the wind in the equatorial plane. 


6.2. The solar wind in time 

Our predictions of the properties of the winds of solar mass 
stars are likely to be more reliable than for lower-mass stars 
given that our wind model is based on the solar wind, and 
that our rotational evolution models are based on better 
observational constraints at solar masses. Although the dis¬ 
cussion in this section is about the time evolution of the 
solar wind, our results apply equally to the winds of other 
solar mass stars. 

In Fig. [lOj we show our predicted mass loss rates 
for the solar wind as a function of time along the 10th 
and 90th percentile rotation tracks. At 100 Myrs, we 
predict that the solar wind mass loss rate would have 
been ~ 7 x 10~ 14 M 0 yr _1 if the Sun had been close 
to the 10th percentile of the rotation distribution, and 
~ 4 x 10 -13 M 0 yr _1 , if the Sun had been close to or above 
the saturation threshold. An interesting result is that the 
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spread in possible mass loss rates for the Sun at 100 Myr 
due to the uncertainties in how rapidly the Sun was rotating 
is only about a factor of six. This is due to the fact that the 
mass loss rates saturate at around 15flg, such that there is 
no difference between the mass loss rates for stars with ro¬ 
tation rates of 15f2g and 100flg. For every track shown in 
Fig. [lO], the total mass losses from the solar wind integrated 
from 100 Myr to 5 Gyr are < 3 x 10 -4 Mg, indicating that 
the winds of solar mass stars are unable to significantly 
influence their masses during their main-sequence lives. 

If we consider just the ages beyond 700 Myr when the 
rotation rates of solar mass stars have converged and follow 
a dependence on age given by Eqn. [7] our model predicts 
that 


M* oc t 


- 0.75 


(23) 


In the right panel of Fig. [TO] we show how our predictions 
for the evolution of the solar wind mass loss rate compare 
with other constraints from the literature. The solid area 
shows the range of wind mass loss rates from our model 
shown in the left panel of Fig. [TOj The solid black lines show 
predictions of the mass loss rate along our 10th and 90th 
percentile ro tation tracks from the model of |Cranmer fe| 
Saar (2011b 1 6 The Cranmer & Saar (2011b) model predicts 
mass loss rates that are a factor of a few higher than the 
mass loss rates that we predict. This is mostly because the 
mass loss rate that we used for the current solar wind is 
approximately a f actor of 2.3 lower than the mass loss rate 
predicted by the Cranmer & Saar (2011b I model, which 
is slightly higher than 3 x 10~ i4 




g yr x As we show in 
Fig. 5 of Paper I, based on spacecraft measurements of the 
solar wind mass flux, this is likely to be too large. The 
black dotted lines in Fig . |T0| show the predictions of the 
Cranmer & Saar (2011b) model scaled down by a factor 
of 2 .3. These p redictio ns mat ch our model very well, with 
the Cranmer & Saar (2011b) model predicting a slightly 
stronger decrease in solar wind mass loss rate with age, 
and therefore a slightly higher mass loss rate at 100 Myr. 

The solid blue line in Fig. [To] shows the solar wind mass 
loss rate scaled backward s in time by assuming M* oc f -1 ’ 23 , 
as predicted by Suzuki et al. (2013). The age dependence is 
slightly steeper than our prediction, but very similar, and is 
in good agre ement with the predictions of the |Cranmer fe| 
Saar (2011bI model. Our m odel could very ea s ily be m ade 
to agree with the m odels of Cranmer & Saar (2011b) and 
Suzuki et al. (20131 if we assumed a slightly weaker scal- 
ing between the dipole field strength and rotation, which is 


certainly plausible. See et al. (2014) assumed the wind tem¬ 


perature is a constant and scaled the wind base density with 
coronal electron densities. Since the wind mass loss rate in 
their model is proportional to the base density, for stars 
with a given radius, their Eqn. 5 implies that M* oc L x . 
Combined wit h the depe nde nce between rotation and X-ray 
emission from Wright et al. (20111 and Eqn. [Tj the model of 


See et al. (2014) implies M* oc t~ 0 62 , which is very similar 
to our determination. 

In Fig. [10] we also compare our predictions to obser¬ 
vational constraints. The downward pointing arrows at the 
top of the figure show upper limits on the mass loss rates 


6 We calcul ate the mass loss rates fro m the model of|Cranmer_fc| 
Saar (201 lb I using the BOREAS code ( Cranmer fe Saar|2011a| ), 
assuming solar metallicity and taking into account the evolution 
of the solar radius and luminosity. 
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Fig. 10. Figures showing the evolution of the solar wind mass loss rate with age on the main-sequence. The left panel shows the 
mass loss rates that we derive from our rotational evolution models, such that Af* oc ST*' 33 for a given stellar radius and mass. The 
lower and upper lines show AT* along the 10th and 90th percentile rotation tracks respectively. The right panel compares these 
estimates with other constraints on the mass loss rates of solar mass stars from the literature. The green shaded region shows the 
results of our w ind model, as shown in the left panel. The solid black lines show the predictions of the model developed by |Cranmer| 


& Saar (2011b), calculated using the BOREAS code, along our 10th and 90th percentile rotation tracks. The dotted black lines 
show the same thing, but with all of the mass loss rates scaled down by a factor of 2.3 in order to make the prediction for the mass 
loss rate of the current solar wind match our value. The solid blue line sh ows t h e current s olar wind mass loss rate extrapolated 


into the past by assuming a time dependence of t 


as predicted by 


Suzuki et al. 


(2013). The three black downward pointing 


arrows show the upper limits on the mass loss rates for three young solar analogues derived from non-detections of radio emission 

by"“ 


Gaidos et al. 


(|2000 l. Th e solid r ed line shows the solar wind extrapolated into the past assuming a time dependence of t 


as suggested by Wood et al. (20051 fr om measurements of as trosph e ric Ly a absorption, and the dashed red line shows the t -1 46 
dependence that we derive in Section 6.2 using the results of Wood et al.Jj 20 051). The red circles show measurements of the mass 
loss rates of several stars using this technique by Wood et al.I £20051, and I Wood et al. (2014). The empty circles correspond to the 


reported mass loss rates normalised by the stellar surface areas, A7*.R7", an d the solid circle s corre spond to Af*J?7 AT* " 3t> , where 
T?* and AT* are in solar units. We have chosen an age of 4 Gyr for e I n d based on King et al. (2010), though significant uncertainty 
exists in the age determination of this star (see Section 7.3 of King et aL|20101. 


for three young solar analogues derived fr om non-detections 
of stellar winds in radio by Gaidos et al. (2000). Clearly all 
other wind mass loss rate predictions shown in Fig. [TO] are 
consistent with these upper limits. More important are the 
measure ments of the mass l oss rates of severa l sola r ana¬ 
logues by Wood et al. (2002), Wood et al. (20051, and Wood 
et al. ( 2014|) . The mass loss rates were derived by measure¬ 
ments oFexcess absorption of the stellar Lya emission line, 
which is thought to be a result of neutral hydrogen walls at 
the edges of the stellar astrospheres. The solid red line in 
Fig. [TO] shows the relation between mass loss rate and age 
oc t -2 ' 33 suggested by 


AT, 


Wood et al. 


(2005). 


The history of the solar wind suggested by Lya mass 
loss rate measurements is in stark contradiction to our pre¬ 
dictions for the mass loss rates. Our results suggest a much 
weaker dependence of mass loss r ate on stellar age than the 
£- 2.33 dependence suggested by Wood et al. (20051, and we 
do not find any break down of the winds at ages younger 
than 700 Myr. Our mass loss rate predictions are a result 
of considering the rotational evolution of solar mass stars, 
and can be calculated based on the dependence of rota¬ 
tional evolution on wind mass loss rates. It is interesting 
to con sider what implications the scenario of’Wood et al.| 
(2005) would have for the rotational evolution of solar mass 


stars. By 700 Myr, the rotation rates of solar mass stars 
have mostly converged onto the slowly rotating track, and 
the stars spin down approximately according to Eqn. [7] 


such that IT* oc t °" 566 , which implies that dVL^/dt oc fI 2 ' 8 . 
Assuming AT* oc t ~ 2 33 implies that AT* oc ST* 12 . Insert¬ 
ing Udip oc ST*' 32 and this dependence into Eqn. 11 gives 
dST*/df oc STf' 5 , which contradicts the Skumanich style spin 
down. 

A possible reason wh y the age depe ndence of the mass 
loss rate derived by Wood et al. (2005) might be overesti¬ 
mated is that they combined stars of different masses to 
derive a single relation between Fx and AT/i? 2 . However, 
our results suggest a further AT* -3 ' 36 dependence for the 
mass loss rate, and so the results could have been biased by 
th e fact that the inte rmediate activity stars in the sample 
of Wood et al. (2005) were mostly K stars, whereas the low 


activity stars were mostly G stars. The four mo st Sun-like 
single stars in the sample of \ 

ffr. If 


(2005) are e Eri, 


I Wo o d et a,l.| 

61 Cyg A, e Ind, and 61 Vir. Including the Sun in this sam¬ 
ple, we find that a power-law fit between F\ and AT/T? 2 
gives 


— oc F 1 ' 33 


(24) 


in agreement with the results of Wood et al. (2005). If we 
instead include the AT“ 3 ' 36 dependence, we find 


RIM. 


M. _ (x pi -18 

-3.36 oc ) 


(25) 
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SLOW WIND 




FAST WIND 



FAST WIND 



Age (Myr) 


Fig. 11. Figure showing the evolution of the solar wind speed (upper row) and density (lower row) at 1 AU on the main-sequence. 
In all panels, the lower and upper lines show the properties along the 10th and 90th percentile rotation tracks respectively. We 
show the evolution of the slow (left column) and fast (right column) components for both Model A (blue) and Model B (red). In 
all panels, the solid lines show the wind properties assuming only thermal pressure driving of the wind and the dashed lines show 
the results for the 90th percentile track at 1 AU in the equatorial plane taking into account magneto-rotational acceleration of the 
wind, as described in Section[6.1| 


which is indeed shallower]^] How ever, co mbining thi s with 
the relation Fx oc t -174 used by Wood et al. {20051 leads 


M oc t~ 1M , 


(26) 


which is much closer to the M oc t 0,75 that we d erive . 
If we instead assume that fl* oc t -0 - 5 (Skumanich 1972) 


and Rx oc fl 2 (Reiners et al. 


2014), we find M oc t 


- 1.18 


to an age dependence of M that is still significantly steeper 
than our prediction. A large part of this contradiction can 
be further resolved if we consider that the Fx oc t -174 
relation is likely t oo str ong. Com binin g S2* oc f -0,566 and 
Rx oc SI 218 from Wright et al. (20111 gives Fx oc i -1 ' 23 , 
which combined with Eqn. [25] gives us a relation for solar 
mass stars of 


In Fig. [lO] we compare Eqn. [26] with the other predictions. 
The previous f -2 ' 33 determination appears to fit e Eri much 
better than the new determination, though it could simply 


be that the age of ~1 Gyr that w e use is overest imated; the 
gyrochronological age derived by Barnes (2007) of 460 Myr 
fits Eqn. [26] much better. 

The suggestion that the relation between M* and age 
breaks down at an age of 700 Myr, with younger stars hav¬ 
ing low mass loss rates, is also difficult to reconcile with the 
observed rotational evolution of stars at these ages. Assum¬ 
ing that young solar analogues have mass loss rates similar 
to that of the current solar wind is similar to assuming in 
our rotational evolution model that a ~ 0 at these ages, 

' The Fx' 18 dependence is highly sensitive to the physical pa¬ 
rameters for t Eri, since it is the only star at intermediate ac¬ 
tivity in t he sample. The radius of 0 .74 Rq that we use was 
derived by Baines & Armstrong ( 20121 by comparing their mea¬ 
sured interferometric radius with the star’s parallax. They used 
their results to derive a mass of 0.82 Mg and age ~1 Gyr from 
the Yonsei-Yale isochrones. 
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where M* ex f2“. The rotation track that we would expect 
for a rapidly rotating star in this scenario is shown as the 
dark green line in the left panel of Fig. [5] In the first few 
hundred million years after the ZAMS, there is almost no 
spin down, in contradiction to the observational constraints. 
For example, observations indicate that a star at the 90th 
percentile of the distribution of rotation periods at 150 Myr 
would have a rotation rate of ~50 Hq, and that it would 
spin down to ^5 Hq by 550 Myr. However, the rotation 
track in Fig. 5] indicates that if we assume the mass loss 
rates are similar to that of the current solar wind, then the 
star would instead only spin down to ~30 CIq. 

In Fig. 0 we show the wind speeds and densities for 
the slow and fast winds calculated from both Model A and 
Model B as a function of age. Although the mass fluxes in 
Model A and Model B are the same, the two models can 
lead to significantly different predictions for the density and 
speed of the wind. In Model A, since we scale the wind tem¬ 
perature with coronal temperatures, rapidly rotating stars 
have much higher wind temperatures than slowly rotating 
stars, and therefore much higher wind speeds. On the other 
hand, the wind temperature for Model B is determined by 
the surface escape velocity, and is therefore independent 
of rotation. For Model B, the solar wind speeds are ap¬ 
proximately constant in time, with small decreases of about 
70 km s -1 between 100 Myr and 5 Gyr. The decrease in the 
wind speed is due to the expansion of the Sun which causes 
the wind base temperature to decrease. For Model A, the 
wind speeds are much larger at young ages, with values 
as high as 2000 km s -1 for the fast wind in the saturated 
regime. In both models, large changes in the proton densi¬ 
ties at 1 AU are seen between the different rotation tracks, 
and between 100 Myr and 1 Gyr. However, the variations in 
the densities in Model A are a factor of a few smaller than 
the variations in the densities in Model B. The differences 
in the densities between Model A and Model B cancel out 
the differences in wind speeds, leading to the same mass 
loss rates in both models. An important difference between 
the two models is therefore the wind ram pressure, which is 
given by pv 2 . Since the wind ram pressure is more sensitive 
to the wind speed than the density, in Model A the ram 
pressures for the young solar wind will be a factor of a few 
higher than in Model B. 

We should however be cautious when interpreting the 
results for the 90th percentile track shown in Fig. 0 For 
the most rapidly rotating stars, magneto-rotational accel¬ 
eration of the wind is likely to dominate over thermal accel¬ 
eration. As we show in Paper I, the rotation rates at which 
solar mass stars transition into the FMR regime depend on 
which model we use to calculate the wind temperature: for 
Model A and Model B respectively, the transition from the 
SMR regime to the FMR regime happens at ~ 60Hq and 
~ 15Hq. Since most solar mass stars never reach such high 
rotation rates, we conclude that the Sun probably did not 
make it into the FMR regime, though there is still a signif¬ 
icant possibility that it did. The dashed lines in the upper 
panel of Fig. 0 show the Michel velocities along the 90th 
percentile rotation track where they become larger than the 
thermal speeds. For a solar mass star at the 90th percentile 
of the rotational distribution at 100 Myr, the terminal wind 
speed in the equatorial plane is likely to be ^2000 km s -1 , 
regardless of the wind temperature. The winds therefore 
have a much lower density in the equatorial plane than we 
would predict from thermal acceleration alone. 
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Fig. 12. Figure showing the distribution of stellar mass loss 
rates at 100 Myr, 500 Myr, 1000 Myr, and 2000 Myr based 
on our scaling law for mass loss rate and the distributions of 
rotation rates shown in Fig. [7] 
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Fig. 13. Figure showing the evolution of the slow and fast wind 
speeds at 1 AU for different stellar masses assuming Model A 
for determining the wind temperature. In this model, we assume 
that the wind temperature can be scaled with the coronal tem¬ 
perature. We calculate the wind speed at 1 AU using Eqn. |19| 
The three panels correspond to different ages, and in each panel, 
the lower and upper lines correspond to the wind speeds for stars 
at the 10th and 90th percentiles of the rotational distributions 
respectively. In the case of the low-mass stars in the upper panel, 
there is only one line because both percentiles lie above the sat¬ 
uration threshold, and therefore have the same wind tempera¬ 
tures. The lower and upper dashed black lines show the Michel 
velocity along the 10th and 90th percentile tracks, which is an 
estimate of the terminal wind speed in the equatorial plane due 
to magneto-rotational acceleration, as discussed in Section [6.1| 
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There is clearly a lot of uncertainty in the properties of 
the solar wind at young ages. This is not just due to the 
lack of theoretical understanding of stellar winds, but also 
because we do not know how fast the Sun was rotating in 
its first few hundred million years. In the following section, 
we extend our discussion to the winds of lower-mass stars. 
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6.3. The winds of stars from 0.4 Mq to 1.1 Mq 

We are now in the position to study the time evolution of 
stellar wind properties on the main-sequence for a range of 
stellar masses. In Fig. [7] we show the time evolution of rota¬ 
tion periods for a distribution of ~1500 stars from 100 Myr 
to 2000 Myr. In Fig. [T2] we show the evolution of the dis¬ 
tribution of mass loss rates for the stars in this sample. 
The mass loss rates for different stellar masses at the 10th, 
50th, and 90th percentiles of the rotational distribution at 
different ages are given in Table [2] 

The mass dependence of the mass loss rate depends 
strongly on the rotation rate. In the unsaturated regime, 
low-mass stars have higher mass loss rates by a factor of a 
few than high-mass stars, at a given rotation rate. Once all 
stars are in the saturated regime, the mass dependence of 
mass loss rate has been reversed, with the low-mass stars 
having the lowest mass loss rates. This is due to the low 
saturation threshold for low-mass stars. Combining Eqn. [4] 
and Eqn. [g] and assuming that I?* oc Mf 8 , gives 

Ms at * 37 M 0 (^J , (27) 

where M sat is the saturation mass loss rate. This means 
that on the main-sequence, the maximum possible mass 
loss rate for low-mass stars is lower than the maximum 
possible mass loss rate for high-mass stars. This also means 
that the saturation mass loss rate for solar mass stars is 
37 Mq, corresponding to 5.2 x 10 -13 M 0 yr _1 (in TableJSJ 
the saturation threshold for solar mass stars at 100 Myr is 
slightly lower than this value because in those calculations 
we take into account the slow increase of the stellar radius 
with age on the main-sequence). 

At 100 Myr, most low-mass stars lie at the saturation 
mass loss rate and only a few have low mass loss rates. The 
distribution at higher masses also has a significant number 
of stars at the saturation threshold. However, for stars with 
masses above 0.7 M 0 , the distribution has a very strong 
cluster of stars with relatively low mass loss rates, corre¬ 
sponding to ~ 10~ 13 M 0 yr _1 , which is not much higher 
than the mass loss rate of the current solar wind. For ex¬ 
ample, for 0.8 M 0 stars, the mass loss rates at the 10th and 
50th percentiles of the rotational distribution are 5M 0 and 
1OM 0 respectively. This is unsurprising given that these 
stars are only rotating a few times faster than the current 
Sun. At 500 Myr and 1000 Myr, the low mass loss rate 
feature is more dominant, with a much smaller number of 
stars lying at the saturation value. By 2000 Myr, no satu¬ 
rated stars remain at any mass, and all of the stars lie on 
the track of low mass loss rates. At this age, the mass loss 
rates of low-mass stars are generally slightly higher than the 
mass loss rates of high-mass stars. This track is influenced 
by two competing effects: on the one hand, low-mass stars 
have higher mass loss rates at a given rotation rate than 
high-mass stars, but on the other hand, low-mass stars ro¬ 
tate slower than high-mass stars. The result is only a weak 
dependence of the mass loss rate on stellar mass at a given 
age. 

In Fig. [13] we show the slow and fast wind speeds at 
1 AU for stars lying on the 10th and 90th percentiles of 
the rotational distributions at 100 Myr, 1000 Myr, and 
5000 Myr, based on Model A. Along the 10th percentile 
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Fig. 14. Figure showing proton densities (left column) and ram pressures (right column) as a function of stellar mass at three 
different stellar ages for the slow wind component. In each panel, blue corresponds to Model A and red corresponds to Model B, 
and the lower and upper lines correspond to the 10th and 90th percentiles of the rotational distributions. 


track, there is almost no dependence on stellar mass. How¬ 
ever, along the 90th percentile track, the wind speeds are 
much higher for high-mass stars due to the higher satura¬ 
tion threshold. The saturation wind speeds for 0.5 M 0 stars 
are ^700 km s” 1 and ~1200 km s -1 for the slow and fast 
wind models respectively. At later ages, due to their slower 
rotation, the wind speeds are lower for low-mass stars. We 
do not show in Fig. [13] the wind speeds calculated from 
Model B since they are approximately constant with age 
and independent of stellar mass. 

The lower and upper dashed lines in each panel of 
Fig. [ 13 ] show the Michel velocities along the 10th and 90th 
percentiles of the rotational distributions. At 100 Myr, 
along the 10th percentile track, the Michel velocity is much 
lower than the thermal wind speeds. On the other hand, 
along the 90th percentile track, the Michel velocity is larger 
than the thermal wind speeds of the slow component, and 
therefore are likely to determine the wind speeds, even in 
Model A. At later ages, the Michel velocity becomes much 


lower than all of the thermal wind speeds, and therefore 
magneto-rotational effects can be ignored. 

In Fig. |14[ we show the time evolution of wind density 
and ram pressure at 1 AU for both Model A and Model B. 
For simplicity, we only show the results for the slow wind 
model, though we find similar trends for the fast wind 
model. As in previous plots, blue corresponds to Model A 
and red corresponds to Model B, and the lower and up¬ 
per lines in all panels correspond to the 10th and 90th 
percentiles of the rotational distributions respectively. The 
influence of magneto-rotational acceleration of the wind 
is shown as the dashed black line in the upper panels of 
Fig. [l4j For both models, the wind densities and ram pres¬ 
sures at 100 Myr are higher than at 5000 Myr due to the 
higher mass fluxes. At 100 Myr, the densities for Model B 
are a factor of a few larger than for Model A due to the 
lower wind speeds. On the other hand, the ram pressures 
calculated from Model A are larger than the ram pressure 
from Model B. This is because the ram pressure, propor- 
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tional to M+v, has a stronger dependence on wind speed 
than the mass flux. 

The evolution of the distributions of wind density and 
ram pressure is interesting for both Model A and Model B. 
At 100 Myr, there is no spread in these properties for low- 
mass stars since they are mostly saturated. This lack of 
spread is of course partly artificial since we assume one uni¬ 
versal relation between mass loss rate and stellar mass and 
rotation, whereas in reality, there is likely to be some scat¬ 
ter about this relation. Going to higher masses, there begins 
to be a spread in the wind densities and ram pressures be¬ 
tween the slowest and fastest rotators and this spread has a 
strong mass dependence. By 1 Gyr, due to the fast conver¬ 
gence of the rotation rates, this spread in wind properties 
has almost disappeared for solar mass stars. On the other 
hand, a spread in the wind properties for low-mass stars 
has developed. By 5 Gyr, since all of the stellar rotation 
rates have converged, there is almost no spread in the wind 
properties at any stellar mass. 

Interestingly, at 5 Gyr, the wind properties do not have 
a strong mass dependence, and due to the slow evolution of 
rotation at later ages, also evolve slowly with age. In studies 
that require only a crude estimate for the properties of an 
individual inactive star’s wind, it would be reasonable as 
a first approximation to simply assume the measured solar 
wind properties, even if the star has a different mass to 
the Sun. This is only true for slowly rotating stars, and 
it is unclear how well this extends to stars with masses 
below 0.4 M 0 . Of course detailed modelling of the winds of 
individual stars should be preferred when possible. 


7. Summary and Discussion 

Stellar rotation and winds are fundamentally linked on the 
main-sequence. Stellar winds remove angular momentum 
from their host stars, causing their rotation rates to de¬ 
crease as they age. At the same time, the decrease in the 
rotation rates of stars likely causes large changes in the 
wind strengths. In Paper I of this series, we develop a stel¬ 
lar wind model based on scaling the solar wind to other 
stars. In this paper, we develop an observationally driven 
rotational evolution model which we couple to our wind 
model. We develop the rotational evolution model for two 
purposes: firstly, we use the model to constrain the wind 
mass loss rates of low-mass main-sequence stars as a func¬ 
tion of stellar mass, radius, and rotation, and secondly, we 
use the results of the model to show how wind properties 
evolve on the main-sequence. 

To constrain our rotational evolution models, we col¬ 
lect measured rotation periods for stars at different ages 
to determine the evolution of the distribution of rotation 
rates on the main-sequence. Probably the major gap in our 
observational understanding of rotational evolution is the 
rotational distribution of low-mass stars (M* < 0.8M 0 ) at 
intermediate ages (~l-2 Gyr). In Fig. [2j we show the obser¬ 
vationally constrained rotational distributions at 150 Myr 
and 550 Myr. Statistics related to these distributions are 
given in Table [T] At 150 Myr, there is a massive spread in 
the rotation periods of stars of all masses, though higher 
mass stars are more concentrated in the slowly rotating 
track. By 550 Myr, the rotation rates of solar mass stars 
have mostly converged, and some convergence has taken 
place for lower mass stars. The reason for this difference 
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is likely to be the lower saturation threshold for mass loss 
rate and magnetic field strength for lower mass stars. 

It is interesting to speculate about where the Sun might 
have been in the distribution of rotation rates at 100 Myr. 
At this age, approximately 70% of solar mass stars are on 
the slowly rotating track. The dominance of the slowly ro¬ 
tating track can be seen from the fact that the 10th and 
50th percentiles are at 4 U 0 and 7 f2 0 respectively. Naively, 
we would therefore expect that the Sun at 100 Myr was a 
slow rotator, with a rotation rate of ~4 days. However, the 
rotational distribution contains a significant rapidly rotat¬ 
ing tail, with the 90th percentile being at 50 fl 0 . There is 
a significant possibility that the Sun was a rapid rotator 
at 100 Myr. This uncertainty is important because these 
scenarios correspond to radically different histories for the 
Sun’s magnetic activity and therefore its levels of high en¬ 
ergy radiation and the strengths of its w inds . It is co mmon 
in the literature to use the scaling laws of Ribas et al. (2005) 
to predict how the Sun’s X-ray, EUV, and FUV radiation 
has evolved in time. This corresponds approximately to the 
assumption that at 100 Myr, the Sun was rotating with a 
period of ~2.7 days (the rotation period of EI< Dra), which 
is at the 75th percentile of the rotational distribution, and 
is therefore a relatively rapid rotator with a relatively high 
level of magnetic activity. If instead, the Sun had a rotation 
period of 5 days, its X-ray, EUV, and FUV fluxes would 
have been significantly lower, likely correspondin g to the 


l evels of emission for n 1 UMa in the sample of Riba s et al. 
(2005), which has an age of ~300 Myr. If the Sun had been 
an extremely rapid rotator at 100 Myr, it would have had 
a level of emission similar to that of EK Dra, given that 
EK Dra is close to the saturation threshold, but this level 
of emission would not have immediately decayed with age. 
Instead, the Sun would have taken several hundred million 
years to spin down to the saturation threshold, and during 
this time, the emission would have remained at the satu¬ 
ration level. Only after a few hundred million years would 
the emission have started to decay. These uncertainties are 
very significant for studies of the evolution of planetary at¬ 
mospheres given the sensitivity of their upper atmospheres 
to the level of high-energy radiation f rom the central star 
( Tian et al.||2008 Lammer et al.pOlO l. 

The basis of our rotational e volutio n mo del is a for¬ 
mula for wind torques derived by Matt et al. (2012) from 
a grid of MHD wind models. The formula relates wind 
torque to the stellar properties, mass loss rates, and mag¬ 
netic field strengths. For the magnetic field strengths, we 
scale the solar magnetic field strength to other stars us- 
ing the scal ing law for the large scale field given by |Vi- 


clotto et al. (2014a). We then specify the mass loss rate by 
M* ex where a ss 1.33 and b w —3.36 are free pa¬ 

rameters in our model that we fit to the observational con¬ 
straints. Another free parameter in our model is the mass 
dependence of the saturation threshold for Af* and Hdi p , 
which we assume is Af 9 ' 3 . We choose this value since it is 
the closest to the weaker mass dependences derived for the 
saturation of X-ray emission. Interestingly, we find that the 
mass loss rate in the saturated regime depends on AT) 3 , 
which is very different to the M~ 3 36 dependence that we 
find in the unsaturated regime. These things can all be put 
together to derive simple scaling laws for the wind torque, 
as we show in Section |4.2| In the unsaturated regime, the 
wind torque is r w oc , and in the saturated regime, 
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the wind torque is r w oc (where the constants of 

proportionality are different for these two scaling laws). 

Given that we have not considered the rotational evolu¬ 
tion of fully convective stars (M+ < 0.35 M 0 ), it is unclear 
how well o ur scaling laws apply in this mass regime. Wood| 
et al. (20011 measured astrospheric Lya absorption for the 
a Centauri system but was unable to achieve a measure¬ 
ment for the much lower mass Proxima Centauri, and sug¬ 
gested that Proxima Centauri has a mass loss rate that is at 
least an order of magnitude lower than the mass loss rate 
in the a Centauri system. Given the parameters of these 
stars, we would expect from our scaling laws that Prox¬ 
ima Centauri has a slightly larger mass loss rate than both 
components of the a Centauri system combined, which is in¬ 


consistent with the results of Wood et al. 


(2001). This could 


be evidence for a fundamental difference'in the properties of 
the winds of fully and partially convective stars. Although 
highly speculative, such an interpretation is consistent with 
the fact that the properties of global photospheric magnetic 
fields appear to differ significantly b etween t hese tw o types 
of stars, both on the main-s equence dM orin et al.|2008 ) and 
on the pre-main-sequence (Gregory et al. 2012), and with 
the significantly different rotational evolution That has been 
suggested for fully convective stars ( Irwin et al.||2011 |. 

A major problem with stellar wind models is the lack 
of observational constraints on wind properties. In the fi¬ 
nal section of Paper I, we discussed promising methods for 
constraining wind properties, such as comparisons between 
astrospheric Lya absorption and radio interferometric mea¬ 
surements of radio emission from the wind, or stellar obser¬ 
vations during planetary transits. The rotational evolution 
of low-mass main-sequence stars can potentially provide the 
most important observational constraints on stellar wind 
properties. The formulae provided in this paper could be 
used as simple estimates for wind torques, which can then 
be used to constrain free parameters, such as the base den¬ 
sity, in more complex MHD wind models, or to check the 
results when constraints on these parameters are already 
available. The scaling laws given above can be used to esti¬ 
mate the wind torques, but they are not likely to be realistic 
for the most rapidly rotating stars. More accurate is to cal¬ 
culate the torque from Eqn. [3] using Eqn. [2] to get t' . Eqn. [5] 
and Eqn. [E] should be us ed as input into Eq n. [il for M t and 
-Bdip (where Eqn. 11 of jWright et al. |2011 should be used 
to calculate the convective turnover times). If the star is a 
rapid rotator, then Eqn. [6] should be used to calculate the 
saturation rotation rate, with c = 2.3. 

In Section [6] we combine the results of our rotational 
evolution model with our wind model developed in Paper I 
to show how wind properties evolve on the main-sequence 
for stellar masses between 0.4 M 0 and 1.1 M 0 . For the so¬ 
lar wind, we find a much weaker depe nde nce of mass loss 
rate on age than predicted by Wood et al. (2005), with the 
highest mass loss rates at young ages that are only a factor 
of 20 above that of the current solar wind. We also find no 
evidence for weak mass loss rates at young ages. The strong 
winds at young ages, and the weak dependence of mass loss 
rate on age at intermediate and late ages, appear necessary 
for solar mass stars to spin down in a way that is consis¬ 
tent with observational constraints on rotational evolution. 
We show that the early solar wind was likely to have been 
denser than the current solar wind, though the magnitude 
of this change depends in most cases on how wind temper¬ 
atures scale to other stars, since temperature determines 


the wind speed. This is not true for the most rapidly ro¬ 
tating stars, given that the wind speeds are determined by 
magneto-rotational acceleration. 

The large spread in rotation rates for young stars leads 
to large spreads in the wind properties, though this spread 
is to some extent decreased by saturation. Interestingly, 
due to the fact that most stars with masses ^0.5 M 0 lie 
above the saturation threshold at 100 Myr, this spread is 
not present at these masses. As stars spin down, the spread 
in wind properties of solar mass stars disappears, but a 
small spread for lower mass stars develops. However, by a 
few Gyr, this spread has completely disappeared due to the 
convergence of rotation rates. 

Our motivation for this study is to provide important 
input into further studies that attempt to improve our un¬ 
derstanding of the complex interactions between stars and 
the evolving atmospheres of potentially habitable planets. 
One interesting question in this area is the way in which 
the current Earth’s atmosphere formed. It is possible that 
the primordial Earth formed with a dense hydrogen atmo- 
sphere captured during the Sun’s classic al T Tauri phase 
(|Hayashi et al.|1979| |Lammer et al.||2014| which must have 
been lost early in the Earth’s evolution. It is easy to imag¬ 
ine how this could have happened if the Sun was born as 
a rapid rotator: in this scenario, the Sun would have re¬ 
mained highly active for a long time, driving significant 
thermal and non-thermal escape from the Earth. However, 
not only was the Earth able to lose this original hydrogen 
atmosphere, but it was able to do so while keeping the thin 
nitrogen domin ated atmosphere that it currently has. |Lichfi| 
enegger et al. (20101 estimated that at ages below ~1 Gyr, 
the combination of a stronger solar wind compressing the 
Earth’s magnetosphere and a higher EUV flux inflating the 
Earth’s atmosphere leads to a significant portion of the at¬ 
mosphere being exposed to the solar wind. They estimated 
that this would have caused so much non-thermal escape of 
atmospheric gas that the current Earth’s atmosphere would 
have been removed in a very short amount of time. A possi- 
ble sol ution to the problem identified by Lichtene gger et al.| 
(20101 is that the Sun was instead a slow rotator in its early 
main-sequence life, and therefore did not have EUV emis¬ 
sion and wind strengths significantly above current levels 
at any time after the ZAMS. 

Although there are potential solutions to these problems 
that do not depend on the early Sun’s level of activity, the 
above considerations show that the rotation history of the 
Sun can be a significant uncertainty in our understanding 
of the formation of the Earth’s atmosphere. The Earth is 
the only planet that we know has formed in such a way 
as to be habitable. The origins and the evolution of stellar 
rotation, and how these connect to winds and high-energy 
radiation, are important factors that need to be considered 
when studying how often this is likely to have happened in 
other stellar systems. 
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